Discrete Biochemical Systems Theory

Almost every biomedical systems analysis requires early decisions regarding the choice of the most suitable representations to be used. De facto the most prevalent choice is a system of ordinary differential equations (ODEs). This framework is very popular because it is flexible and fairly easy to use. It is also supported by an enormous array of stand-alone programs for analysis, including many distinct numerical solvers that are implemented in the main programming languages. Having selected ODEs, the modeler must then choose a mathematical format for the equations. This selection is not trivial as nearly unlimited options exist and there is seldom objective guidance. The typical choices include ad hoc representations, default models like mass-action or Lotka-Volterra equations, and generic approximations. Within the realm of approximations, linear models are typically successful for analyses of engineered systems, but they are not as appropriate for biomedical phenomena, which often display nonlinear features such as saturation, threshold effects or limit cycle oscillations, and possibly even chaos. Power-law approximations are simple but overcome these limitations. They are the key ingredient of Biochemical Systems Theory (BST), which uses ODEs exclusively containing power-law representations for all processes within a model. BST models cover a vast repertoire of nonlinear responses and, at the same time, have structural properties that are advantageous for a wide range of analyses. Nonetheless, as all ODE models, the BST approach has limitations. In particular, it is not always straightforward to account for genuine discreteness, time delays, and stochastic processes. As a new option, we therefore propose here an alternative to BST in the form of discrete Biochemical Systems Theory (dBST). dBST models have the same generality and practicality as their BST-ODE counterparts, but they are readily implemented even in situations where ODEs struggle. As a case study, we illustrate dBST applied to the dynamics of the aryl hydrocarbon receptor (AhR), a signal transduction system that simultaneously involves time delays and stochasticity.

Introduction

Arguably the greatest challenge of systems modeling in the biomedical sciences is the choice of optimal process representations. Often the true magnitude of this challenge is ignored and the modeler either constructs an ad hoc model or chooses a default, such as a Lotka-Volterra system for describing the interactions among competing populations (Volterra, 1926; Lotka, 1956; May, 1973) or a mass action formulation or some variation of the Michaelis-Menten rate law for enzyme catalyzed processes (Michaelis and Menten, 1913; Voit et al., 2015). These default representations may be further extended or refined with the inclusion of environmental variables in a population model (Stein et al., 2013; Dam et al., 2020) or the inclusion of modulating effects, such as the regulation of a biochemical reaction through competitive or allosteric inhibition (Cornish-Bowden, 2012). Because there are no iron-clad rules for choosing a model, researchers often arrive at rather different formulations even for the same phenomenon. An illustrative example is the phosphofructokinase reaction in glycolysis, for which numerous rate functions of drastically different complexity have been proposed (Voit, 2017a). The choice of optimal representations becomes even more challenging at the intersections of typical biological domains, such as the combination of genetics, metabolism, and organismal physiology, because the default models of the various subdisciplines are different, thereby creating the need of multiscale models that operate at different temporal, spatial and organizational scales.

One could argue that biological processes must obey the laws of physics and that, therefore, optimal—or at least adequate—representations are prescribed. While this is true in a fundamental sense, most biological processes are so convoluted that exact physical representations of all contributing aspects become infeasible (Voit, 2008). As an example, consider the generation of two daughter cells from a bacterial mother cell. At a high level, one bacterium becomes two, two become four, and so on, and it is easy to formulate an exponential function that describes the progression well. However, if it is necessary to account for more details, for instance, in order to understand a mutant with aberrant behavior, it becomes clear that the cell division process is immensely complicated (Schafer, 1998; Carlton et al., 2020). It is multifaceted and involves so many different aspects at the molecular level that it is hardly possible to formulate the governing processes, proceeding in time and space, with elementary functions that are directly derived from the first principles of physics.

A second aspect of the challenge of biomedical systems modeling is the fact that it is usually difficult to capture the dynamics of a molecular or cellular component directly. Even the simple Michaelis-Menten rate law of enzyme kinetics does not prescribe the changing concentration of a substrate or product as the reaction progresses, but expresses the speed of the reaction as a function of the substrate concentration. By contrast, it is often feasible to characterize all influences that lead to an increase or decrease in a system component over time (Voit, 2020). Indeed, the literature contains uncounted articles about “the effect of … on …,” which explicitly or implicitly describe how a target variable changes in response to some input. Thus, this view focuses on the change in a component, rather than the state of this component, and this change is driven by the totality of all contributing factors. A natural mathematical formulation of this situation is a system of ordinary differential equations (ODEs) which, after all, equate the instantaneous change in a variable to all processes affecting this variable. Consequently, the biomathematical literature contains an enormous body of work using ODEs to analyze biological systems [for introductory texts, see (Keshet, 2005; Klipp et al., 2016; Voit, 2017b)]. Even so, it must be kept in mind that ODEs are approximations of natural processes, which are often genuinely discrete (see Supplementary Data S1, S2).

While ODEs have become the standard modeling default, the conundrum of determining the best possible model structure persists. Two generic solutions are 1) the use of ad hoc representations that are often chosen simply for convenience and match the natural processes sufficiently well and 2) suitable, unbiased approximations. Among the latter, linear systems are most straightforward but are often at odds with the genuine nonlinearities of biomedical systems. A prominent alternative is Biochemical Systems Theory (BST) (Savageau, 1976; Voit, 2000; Torres and Voit, 2002; Voit, 2013), which uses power-law representations for all processes, thereby creating highly structured nonlinear models in immutable, predefined formats (Supplementary Data S1).

Independent of what representations are chosen to design ODE models, the ODE format in itself faces a number of challenges. Of particular prominence among these are time delays and stochastic effects (Supplementary Data S2). Sometimes, these can be addressed with sophisticated numerical ODE solvers, but the formulation and implementation can quickly become convoluted and often requires intimate knowledge of the inner workings of these solvers.

An illustrative example for the crucial role of delays is a situation that arose when we analyzed the dynamics of anemia during malaria, a disease that is caused by Plasmodium parasites that invade red blood cells and eventually cause them to burst. Red blood cells furthermore disappear in large numbers due to a so-called bystander effect, in which many non-infected red blood cells perish for reasons that are not well understood. One difficult challenge that arose during our modeling attempts was the fact that red blood cells naturally have a narrowly determined life span with rather small variation; in humans, it is about 115 days ± 15% (Franco, 2012). The modeling challenge becomes apparent in the assessment of how many cells are expected to disappear at a given time point during the infection (Fonseca and Voit, 2015; Fonseca et al., 2016). Some disappear due to the infection or the bystander effect, but many are removed by the spleen because they have reached the end of their natural life. To account for the latter aspect, one needs to know the age of each cell at any given point in time. However, ODEs do not recount the ages of individual cells. Thus, the disappearance of cells from the blood stream must be based on averages, which are adequate under steady-state conditions, but not for dynamic changes caused by the growing parasite population. Even the use of delay differential equations (DDEs) is inconvenient in this case, whereas a discrete, recursive modeling approach is straightforward (Fonseca and Voit, 2015; Fonseca et al., 2016). Other pertinent examples of delays and stochasticity are presented in Supplementary Data S2.

This article proposes an alternative to BST models that facilitates the modeling of genuine discreteness, delays, small numbers of components, stochastic events and combinations of these complicating factors. This alternative consists of a discrete, recursive version of BST, here dubbed “dBST,” which is straightforwardly constructed and implemented.

Results

For the practicing computational modeler in the biosciences, a partial solution to the drawbacks of ODEs can be the use of systems of discrete-time, recursive equations, where the changes in variables are represented on the basis of power-law functions, as is the case in BST. This replacement of ODEs with recursive equations raises the immediate question whether any genuine features of ODE models are lost. The answer can be approached in two ways. First, it is rather evident that the recursive equations converge to the ODEs if the step size decreases to 0 in the limit. In fact, computer algorithms for solving ODEs use small discrete step sizes. Second, one may test whether representative nonlinear phenomena that are typically represented with ODEs, such as saturation, limit cycles, and deterministic chaos, can also be represented through recursive equations with a reasonable step size. Supplementary Data S3 discusses mathematical similarities between BST and dBST systems. Here, we focus on the response repertoire of dBST systems and their features. We also present a case study illustrating the de novo design of a dBST system that simultaneously accounts for both, delays and stochasticity.

Response Repertoire of Discrete Biochemical Systems Theory Models

Simple Introductory Example of a Linear Pathway

The 2-variable BST system

${\dot{X}}_1=2X_0-r{X}_1^{0.8}$

${\dot{X}}_2= r{X}_1^{0.8}-2.5X_2^{0.5}(1)$

X₀ = 1

X₁(t₀) = 1.18

X₂(t₀) = 0.64

represents a simple linear pathway with constant input X₀ = 1 and rate r = 1.75 for the conversion of X₁ into X₂. The pathway is shown in Figure 1:

FIGURE 1

FIGURE 1. Simple linear pathway with constant input X₀.

As an illustration, the system is initiated very close to its steady state (1.181653, 0.64). At t = 5, the input X₀ is persistently increased by 20%. Solving the equations shows that the system responds to the changed input by approaching a new steady state (1.484, 0.922) (Figure 2).

FIGURE 2

FIGURE 2. Comparison between the results of corresponding 2-variable BST (lines) and dBST (dots) models (Eqs 2, 10).

The corresponding dBST system in standard notation (Supplementary Data S3) is

$X_{1,q\cdot \vartheta }=X_{1,\left(q-1\right)\cdot \vartheta }+\vartheta \left[2X_0-r{X}_{1,\left(q-1\right)\cdot \vartheta }^{0.8}\right]$

$X_{2,q\cdot \vartheta }=X_{2,\left(q-1\right)\cdot \vartheta }+ \vartheta \left[r{X}_{1,\left(q-1\right)\cdot \vartheta }^{0.8}-2.5X_{2,\left(q-1\right)\cdot \vartheta }^{0.5}\right](\mathrm{2a})$

To simplify this notation for easier reading, we rename ${\tilde{X}}_i=X_{i,q\cdot \vartheta }$ and $X_i=X_{i,\left(q-1\right)\cdot \vartheta }$, for i = 1, 2, which simplifies the appearance of Eq. 2a to

${\tilde{X}}_1=X_1+\vartheta \left[2X_0-r{X}_1^{0.8}\right]$

${\tilde{X}}_2=X_2+\vartheta \left[ r{X}_1^{0.8}-2.5X_2^{0.5}\right](\mathrm{2b})$

Choosing as step size $\vartheta$ = 0.5 reveals output that is quite similar to that of the BST system, although one notes that the responses, especially of X₂, are slightly different immediately following the switch in input (t = 5). Importantly, both formulations exhibit essentially the same dynamics and approach exactly the same steady state (Figure 2). Other step sizes yield similar results.

Limit Cycles

Limit cycles are representations of oscillations that are stable in a sense that, when perturbed by external influences, return to the original frequency and amplitude. Limit cycles are ubiquitous in biology (Keshet, 2005).

It has been proposed that many disease patterns can be seen mathematically as shifts from physiological to pathological limit cycles (Claude, 1995).

Like BST systems, dBST systems can capture the dynamics of stable limit cycles (Lewis and Voit, 1991; Yin and Voit, 2008). An example is the stable oscillator

${\dot{X}}_1=0.011 \left[X_1^2{X}_2^3-X_1{X}_2\right]$

${\dot{X}}_2=0.01\left[X_1^{-1}{X}_2^4-X_1^3{X}_2^5\right](3)$

which in standard dBST format reads

$X_{1,q\cdot \vartheta }=X_{1,\left(q-1\right)\cdot \vartheta } + \vartheta \cdot 0.011 \cdot \left[X_{1,\left(q-1\right)\cdot \vartheta }^2 X_{2,\left(q-1\right)\cdot \vartheta }^3-X_{1,\left(q-1\right)\cdot \vartheta } X_{2,\left(q-1\right)\cdot \vartheta }\right]$

$X_{2, q\cdot \vartheta }=X_{2,\left(q-1\right)\cdot \vartheta } + \vartheta \cdot 0.01 \cdot \left[X_{1,\left(q-1\right)\cdot \vartheta }^{-1}{X}_{2,\left(q-1\right)\cdot \vartheta }^4-X_{1,\left(q-1\right)\cdot \vartheta }^3{X}_{2,\left(q-1\right)\cdot \vartheta }^5\right](\mathrm{4a})$

To the human eye, this format may look rather unwieldy, but it is easily implemented into computer code. Furthermore, using the simplified notation introduced in Eq. 2, we obtain

${\tilde{X}}_1=X_1+\vartheta \cdot 0.011 \cdot \left[X_1^2{X}_2^3-X_1{X}_2\right]$

${\tilde{X}}_2=X_2+\vartheta \cdot 0.01 \cdot \left[X_1^{-1}{X}_2^4-X_1^3{X}_2^5\right](\mathrm{4b})$

Solving the system with step size $\vartheta$ = 0.1 confirms that the system indeed has a stable limit cycle. Namely, initial conditions inside the limit cycle, like (X_1,0, X_2,0) = (1, 1.3), lead to increasing oscillations, while conditions outside, such as (X_1,0, X_2,0) = (1.1, 1.8), generate damped oscillations; starting essentially on the limit cycle, e.g., (X_1,0, X_2,0) = (0.7249193, 0.8685822) demonstrates constant amplitudes. In the phase plane, the corresponding plots are outward and inward spirals, as well as the stable orbit (Figure 3). If we use the step size $\vartheta$ = 1, the system still displays a limit cycle of similar shape, but with larger amplitudes (not shown). Much larger step sizes eventually become too coarse and destroy the features of the limit cycle.

FIGURE 3

FIGURE 3. The dBST system in Eq. 4 models a stable limit cycle, as confirmed by simulations starting inside (A), outside (B) and essentially on the limit cycle (C). Panel (D) displays a typical phase-plane plot with superimposed oscillations spiraling out (dark green) or in (light green) toward the limit cycle, as well as starting very close to the limit cycle itself (cyan); the initial locations are indicated by circles. The trajectories appear to be smooth because the step size is rather small. The ODE model produces essentially the same solutions, even though the maximal amplitudes are slightly different.

Deterministic Chaos

Discrete BST systems are also rich enough to permit deterministic chaos. While a rigorous proof is difficult, an example is a discrete system gleaned from the well-known Lorenz oscillator (Lorenz, 1963), which mathematician and meteorologist Edward Lorenz developed as a simplified representation of atmospheric convection, which had been modeled previously as a fluid layer for which the temperatures at the top and the bottom were kept constant at different values (Saltzman, 1962). The ODE format of this system reads:

${\dot{X}}_1=0.2\cdot X_2-0.2\cdot X_1$

${\dot{X}}_2=0.6\cdot X_1-0.02\cdot X_2-0.02\cdot X_1\cdot X_3(5)$

${\dot{X}}_3=0.02\cdot X_1\cdot X_2-0.05\cdot X_3$

The reformulation into recursive equations is straightforward, and the dBST in simplified notation (see Eq. 2) reads

${\tilde{X}}_1=X_1+\vartheta \cdot 0.2\cdot \left[X_2-X_1\right]$

${\tilde{X}}_2=X_2+\vartheta \cdot \left[0.6\cdot X_1-0.02\cdot X_2-0.02\cdot X_1\cdot X_3\right](6)$

${\tilde{X}}_3=X_3+\vartheta \cdot \left[0.02\cdot X_1\cdot X_2-0.05\cdot X_3\right]$

Solving the system with step size $\vartheta$ = 1 and initial conditions (10, 10, 30) reveals dynamics similar to that of the chaotic Lorenz equations in ODE format, both in the time domain and in phase plane (Figure 4), although the numerical details of the results are different, which is not surprising, as chaotic systems are extremely sensitive to all numerical settings, such as parameter values, initial conditions, and the step size for solving the ODEs or the discrete equations. One also notes that the maximum amplitudes of the BST system are somewhat higher than in the ODE model.

FIGURE 4

FIGURE 4. The dBST system in Eq. 6 captures deterministic chaos, similar to the ODE system proposed by Lorenz. The top panel shows results in the time domain, in comparison to the corresponding ODE system (thin grey lines). The BST and dBST systems diverge quickly (top panel), which is a genuine feature of chaotic systems. The bottom panel displays phase-plane plots of the two models, showing the discrete nature of the system in the form of connected straight lines. The initial locations are indicated with circles. Note that the maximal amplitudes of the dBST system are larger.

Typical Simulations in dBST That Are More Intuitive Than in ODE Models

The simulations described in this section are straightforwardly implemented in dBST, and while it is possible to implement some of them in ODEs, such an implementation is sometimes cumbersome or difficult to intuit. Indeed, one may have to be creative if some of these issues are to be included into ODE solutions and understand the inner workings of the numerical solution algorithms. For example, the widely used Runge-Kutta method averages the slope for each solution step, and additional statements, such as if-conditions, can influence this average or cannot be taken into consideration, depending on how the solver was coded. Furthermore, using numerical solvers with variable step size requires care so that the choice of the optimal step size is not affected.

In the following, we focus on different types of stochastic events and the dependence of the system dynamics on thresholds for dependent variables. Details regarding delays are discussed in Supplementary Data S2 and in the later Case Study Aryl Hydrocarbon Receptor Signal Transduction.

Stochastic Variations in Rates

Returning to the introductory at the beginning of the Results section, it is easy in dBST to replace the constant rate r of the process converting X₁ into X₂ with a rate that stochastically varies within a range of, say, ± 10% of the nominal value in the example. For this illustration, we randomly sampled a rate from this range at every iteration. Two solutions are shown in Figure 5. Variations on this theme are also readily implemented. For instance, it is possible to sample a new rate less frequently than at every step.

FIGURE 5

FIGURE 5. Comparison of results from the deterministic and two instances (A and B) of stochastic models (ϑ = 0.5) of the simple pathway in Figure 2. The two dBST simulations were obtained with a stochastically varying rate r in Eq. 1 (dots), starting from different seeds, while the solid lines are the results of Eq. 2 with constant rate r.

Events Where System Variables Affect External Events

Suppose a signaling system responds stochastically to an environmental trigger, which is a ubiquitous situation in biological systems, especially if the concept of an “environment” includes the biophysical surroundings of cells. We develop this example in several steps, because some cases are easily addressed with ODEs, whereas others are not. In the simplest case of a stochastic environmental input it is possible to include if-statements into a numerical solver, such as the deSolve R package (Soetaert et al., 2010), which was designed to solve various initial-value problems, differential algebraic equations and partial differential equations. However, if one or more of the system variables influence the random variable, the situation is much more complicated, as the random variable must be adressed inside the solver, which is difficult for a ODE solver but straightforward in the discrete case.

Suppose at first that the environmental trigger is present or absent for stochastically long time periods that begin at random time points and whose magnitude affects the response of the signaling system. As a specific example, consider the lac operon of the bacterium E. coli, where external lactose triggers changes in gene expression (Lewis, 2005). Savageau (2001) proposed a model of the system in the form of the diagram in Figure 6 and represented it with S-system equations. In this model, X₁ is the concentration of mRNA of the lac operon, X₂ is the concentration of the enzyme β-galactosidase, which catalyzes the conversion of lactose into galactose and glucose, and X₃ and X₄ are the intracellular and extracellular concentrations of lactose, respectively. X₄ is considered an independent variable and therefore that does not require its own differential equation.

FIGURE 6

FIGURE 6. Diagram of the lac operon and corresponding equations, as proposed by Savageau (2001).

Savageau discussed the kinetic orders (g and h parameters) but did not provide specific parameter values for them or for the rate constants (α and β parameters). We use this example for a series of demonstrations, specifying the parameter values as shown in Eq. 7.

${\dot{X}}_1=0.2X_3^2-0.1X_1$

${\dot{X}}_2=0.5X_1-X_2(7)$

${\dot{X}}_3=0.1X_4-0.1X_2{X}_3$

For the first demonstration, suppose that the stimulus, external lactose (X₄), is available in irregular time periods and concentrations that vary randomly within reasonable ranges. If the switch points and magnitudes of X₄ are known beforehand, the simulation of the ODE system is straightforward. For instance, if the system starts at its steady state (0.43, 0.22, 0.46) with X₄ = 0.1 and the stimulus changes according to Table 1, one directly obtains the output in Figure 7A.

TABLE 1

TABLE 1. A priori known schedule of switches in stimulus in the lac operon model.

FIGURE 7

FIGURE 7. Comparison of responses of BST (faint grey lines) and dBST models (dots) to perturbations in external lactose concentraion (solid lavendar line). (A) Timing and magnitudes are known a priori (Table 1). (B) Signaling events occur in a stochastic manner and signal magnitudes are random within a given range (see Text for details). (C,D) Two sets of responses, where switches in the equation of X₁ depend not only on the stochastic input X₄, but also on the value of X₃. (E) The value of X₃ is used to generate a success probability for a Bernoulli random variable. Specifically, if the Bernoulli process returns 1, the new value for X₄ is given by a truncated normal with mean 0.5 and a standard deviation equal to the value of X₁.

If the timing and magnitudes are stochastic, a simulation with a standard ODE solver can be cumbersome as one needs to know how to evaluate functions of time inside the algorithm. Nonetheless, this situation can still be addressed, for instance, with deSolve or in Matlab. This case is again straightforwardly implemented in a dBST system. One set of results is shown in Figure 7B for event times sampled from an exponential distribution with a rate of 1/60. The magnitude of the signal at every event was sampled from a normal distribution with mean μ = 0.5 and standard deviation σ = 0.25.

As the next phase of the example, we analyze the situation discussed by in Savageau (2001), where the format of the first ODE in Figure 6 depends on the current value of X₃. Specifically, the author defined

${\dot{X}}_1=\{\begin{array}{cc}{\alpha }_{1L}-{\beta }_1{X}_1& if{X}_3<X_{3L}\\ {\alpha }_1{X}_3^{g_{13}}-{\beta }_1{X}_1& if{X}_{3L}\le X_3\le X_{3H}\\ {\alpha }_{1H}-{\beta }_1{X}_1& if{X}_3>X_{3H}\end{array}(8)$

where X_3L and X_3H are threshold values and the corresponding “low” and “high” rate constants α_3L and α _3H are different. While it is possible to address this task by embedding if-conditions into an ODE solver, these situations of thresholds are much more easily called up in dBST: The If-statements are directly implemented in the recursive step for variable X₁. As a demonstration, suppose again that the external trigger changes in unpredictable patterns, which we assume to be random in terms of timing and magnitude, as before, and that the thresholds are in effect. Two typical results are shown in Figures 7C,D.

As the most complicated variation of the example, let us now suppose that both the magnitude and frequency of the stochastic events depend on the state of the system. For instance, the amount of external lactose X₄ to be imported into the cell could stochastically depend on the current mRNA prevalence X₁ and also the internal lactose concentration X₃, which together could have an effect on the characteristics of the import transporters. This situation cannot easily be addressed with an ODE solver, if at all, as one can no longer first calculate the current level of X₄ and then present it to the solver as a time-dependent function. Instead, such a situation mandates that the stochastic variables be evaluated inside the solver, a process that can interfere with the ODE solution. By contrast, the discrete version is easily implemented. An example is presented in Figure 7E.

Parameter Estimation

Much has been written about the estimation of parameter values of ODE systems from time series data [e.g., see reviews (Gennemark and Wedelin, 2007; Chou and Voit, 2009; Gennemark and Wedelin, 2009; Gábor and Banga, 2015)]. In order to make typical gradient methods and evolutionary algorithms more efficient, it has been suggested to smooth the time series data, use points from the smoothed trends, and estimate slopes of the time trends corresponding to the chosen points (Varah, 1982; Voit and Savageau, 1982; Voit and Almeida, 2004). Substituting these quantities in the ODE system converts the task of estimating parameter values for ODEs into an estimation from algebraic equations. Thus, for estimation purposes, the ODE for each variable X_i,

${\dot{X}}_i=F_i\left(X_1, \dots , X_n\right)\mathit{i = 1, ..., }n(9)$

is converted into a system of K algebraic equations of the format

$S_i\left(t_1\right)=F_i\left(X_1\left(t_1\right), \dots ,X_n\left(t_1\right)\right)$

$S_i\left(t_2\right)=F_i\left(X_1\left(t_2\right),\dots ,X_n\left(t_2\right)\right)(10)$

$S_i\left(t_K\right)=F_i\left(X_1\left(t_K\right), \dots ,X_n\left(t_K\right)\right)$

Here, each equation corresponds to one chosen time point. The X values are either the raw data or the corresponding data from the smoothed time trend, while S indicates the corresponding slope of the time course.

This method of using estimated values and slopes tends to be computationally much more efficient than parameter inferences directly from ODEs, for instance with a gradient method or an evolutionary algorithm (Voit and Almeida, 2004). One drawback is that the estimation of slopes exacerbates noise in the data (Knowles and Renka, 2014; Voit, 2017b). To some degree, this problem is alleviated by smoothing the data appropriately.

For the discrete system, no slopes need to be estimated as the difference to be used instead, $X_{i,q\cdot \vartheta }-X_{i,\left(q-1\right)\cdot \vartheta }$, is directly obtained from the data. Thus, given measurements for all X_i at different time points, and possibly a smoothing step, the estimation of the parameters of a dBST system is straightforward.

As an example, consider the branched pathway in Figure 8, for which we pretend to have experimental measurements that had been smoothed, for instance, with a spline. For the illustration, we actually created synthetic “data” from a GMA (BST) model in ODE format and did not worry about noise, in order to assess most clearly to what degree BST and dBST models correspond and reflect the synthetic data. Analogously to the BST model, the format of the dBST equations is dictated directly by the flow structure and regulation of the pathway. In the simplified notation of Eq. 2b, the dBST equations take the form

${\tilde{X}}_1=X_1+\vartheta \cdot \left[a_1\cdot X_3^{g_1}-b_1\cdot X_1^{h_{11}}-c_1\cdot X_1^{h_{12}}\cdot X_3^{h_{13}}\right]$

${\tilde{X}}_2=X_2+\vartheta \cdot \left[b_1\cdot X_1^{h_{11}}-b_2\cdot X_2^{h_2}\right]$

${\tilde{X}}_3=X_3+\vartheta \cdot \left[b_2\cdot X_2^{h_2}-b_3\cdot X_3^{h_3}\right](11)$

${\tilde{X}}_4=X_4+\vartheta \cdot \left[c_1\cdot X_1^{h_{12}}\cdot X_3^{h_{13}}-c_2\cdot X_4^{h_4}\right]$

FIGURE 8

FIGURE 8. Simple branched feedback with dual regulation by X₃.

and we suppose that the values of the parameters are unknown. We analyze datasets with different densities of observation time points. For each parameter optimization, we use the optim function in R (R Core Team, 2018), which is based on the Nelder-Mead method (Nelder and Mead, 1965). Multiple sequential optimizations were performed for each example and the process was stopped when the difference of consecutive errors was less than 10^–3.

For the first illustration, we suppose that data had been obtained in intervals of τ =1, that is, for t = 0, 1, 2, ..., 60, and define $\vartheta$ = 1. The estimation result, shown in Figure 9A, captures the data well. The associated residual error, divided by the total number of data in the four time courses (4 n) is SSE/4 n = 0.1007525/(4 * 61) = 4.13 × 10^–4. Fitting the same data, but with step size $\vartheta$ = 0.5 (results not shown), the fitting error is roughly halved, with SSE/4 n = 0.0443176/(4 * 61) = 1.81 × 10^–4.

FIGURE 9

FIGURE 9. Data fits with dBST [τ = 1 (A) and τ = 3 (B)] and BST (C) models.

As a second illustration, we assume that the data are much sparser (τ =3), that is, with measurements obtained at time points t = 0, 3, 6, 9, … , 60; we again define $\vartheta$ = 1. The result is shown in Figure 9B. The model still fits the data well, with a residual error, divided by the total number 4 n of data, of SSE/4 n = 0.03697914/(4*21) = 4.40 × 10^–4, which is slightly, but not substantially higher than for the denser dataset. The estimated parameter values are slightly different from those obtained for the denser dataset, which is not surprising. However, it is interesting that the sparsity of the data hardly seems to affect the estimation.

Quasi as a baseline for comparison, we also fit the synthetic data with ODE equations in GMA format. They also recapture the data well (Figure 9C), even though the estimated parameter values are not identical to those used to create the data (Table 2), indicating some numerical redundancy among the parameters. The residual error, divided by the number of data is SSE/4 n = 0.03966285/(4*61) = 1.62 × 10^–4, which is again in the same range as for the discrete model. The parameter values are slightly different from those estimated with the dBST model. This result is to be expected because the meaning of each multiplicative parameter is, strictly speaking, not identical for BST and dBST models, as the former represent instantaneous rates and the latter stepwise changes.

TABLE 2

TABLE 2. Parameter estimates obtained for different settings of the dBST model in Eq. 11 and the corresponding BST model.

Case Study: Aryl-Hydrocarbon Receptor Signal Transduction

The aryl-hydrocarbon receptor (AhR) is a highly conserved sensor for specific cues during development and normal physiology (Stockinger et al., 2014; Brinkmann et al., 2019; Zhu et al., 2019), as well as for external, xenobiotic compounds (Stevens et al., 2009; Simon et al., 2015) or danger signals derived from the invasion of parasites, which are mediated through compounds like the tryptophan-derivative kynurenine (Julliard et al., 2014; Gupta et al., 2022). In response to such signals, the AhR signal transduction system triggers the upregulation of a host of genes, most prominently those coding for cytochrome P450 enzymes that metabolize toxicants.

The generic functionality of the AhR-system is depicted in Figure 10. Once a ligand (L) binds to AhR, the activated AhR forms a complex with the AhR nuclear translocator ARNT. This complex translocates to the nucleus, where it serves as a transcription factor that binds to the xenobiotic response element XRE—or a non-canonical XRE analog—within the promotor regions of numerous inducible target genes (Huang and Elferink, 2012). The AhR repressor AhRR competes with AhR for ARNT (Evans et al., 2008). Intriguingly, the gene coding for AhRR is itself under the control of the AhR-ARNT transcription factor, thereby creating a negative feedback loop that eventually stops the expression of AhR-ARNT controlled genes (Zudaire et al., 2008). As one might expect, reality is more complicated, for instance, due to compounds like the hypoxia inducible factor-1α (HIF1α) that compete with AhR and AhRR for ARNT (Spence et al., 1970) and to several cofactors modulating the process (Simon et al., 2015), but the AhR-ARNT-AhRR system by itself contains enough interesting complexity for the present illustration.

FIGURE 10

FIGURE 10. Diagram of the AhR signal transduction system. AhR is activated by a ligand L and binds to the nuclear translocator ARNT. The complex serves as a transcription factor of genes whose promoter regions contain the xenobiotic response element XRE. These genes code for a variety of target proteins (TP) including, notably, the AhR repressor AhRR. AhRR competes with AhR for ARNT, and the complex inhibits gene expression. Transcription and translation incur delays (τ₁, τ₂) and are stochastic in nature (σ₁, ..., σ₄). L is also considered to be stochastic.

One issue in setting up an ODE model is the substantial time delay between transcription factor binding, the actual availability of AhRR, and the resulting repression of target gene expression (Koussounadis et al., 2015). In yeast and mouse, this type of delay was found to be at the order of 3–6 h (Fournier et al., 2010) (Liu et al., 2016). A delay of this magnitude is crucial in the AhR system, as it noticeably delays the inhibitory effect of AHRR on target gene expression.

A second issue is the fact that transcription and translation are known to be stochastic processes (Raj and van Oudenaarden, 2008). In fact, at least in some cases, activation of a promotor causes the production of proteins to occur in short bursts and yields variable protein amounts that occur at random time intervals (McAdams and Arkin, 1997). Delays and stochasticity of course are not mutually exclusive but occur at the same time (Gedeon and Bokes, 2012). A model of this stochasticity for the case of AhR, using a discrete model based on the Gillespie algorithm, was presented by Simon et al. (2015). However, it did not explicitly account for the time delays between AhR binding and target protein expression and the role of AhRR as repressor.

Taken delays and stochasticity into account, we can formulate a dBST model that allows us to test the effects of delay and stochasticity. In mass-action and power-law format, and with simplified notation (see Eq. 2b),such a model has the following format:

The variable names are defined in Figure 10. For this illustration, we choose reasonable rate constants as shown in Table 3 and set the inhibition parameter as g = −4; to avoid numerical issues for X₆ = 0, we define the inhibition as ${\left(X_6+1\right)}^g$. The simulations start at the steady state in the absence of ligand (L = 0), which is (X₁₀, ..., X₉₀) = (10, 0, 10, 0, 0, 0, 0, 0, 0). The shaded terms are subject to delay, stochasticity, or both (see Figure 10). Time is roughly in hours. The step size was taken as $\vartheta$ = 0.1.

TABLE 3

TABLE 3. Parameter values for the AhR signal transduction system.

We show the result of three scenarios. In the first, the time delays and any stochasticity are simply ignored (Figure 11A). The second scenario accounts for two types of delays (Figure 11B), one for transcription and one for translation and activation of protein, and the third incorporates both, delays and noise (Figure 11C). For simplicity, we assume the same delay for the transcription (τ₁ = 3 h) and translation and activation (τ₂ = 4 h) of AhRR and a representative target protein (TP), even though these delays are in reality protein-specific (Koussounadis et al., 2015). Also for simplicity, we assume the same stochastic structure for σ₁ and σ₂ and for σ₃ and σ₄ (see Figure 10). Specifically, these stochastic events are modeled with values from the normal distribution N (1, 0.1), which are multiplied to the affected fluxes. We also added stochasticity to the ligand availability; it did not have much effect but shows up in the dynamics of X₁, ..., X₄.

FIGURE 11

FIGURE 11. Simulation results of three scenarios. (A) Time delays and stochasticity are simply ignored. (B) Delays for transcription and for translation and activation of protein are taken into account. Note that the dynamics of X₇ and X₈ is the same. (C) Ligand availability, transcription and translation are considered stochastic. See Text for further details.

Both simulations start at the steady state without ligand (L = 0). At time t = 2, the ligand concentration is set to 2, and at t = 12 it is returned to 0. The result of the first scenario simulation (Figure 11A) reveals that the production of target protein (TP) very briefly peaks, but that not much TP is produced, due to the immediate onset of inhibition by AhRR. By constrast, accounting for time delays yields a dramatically different dynamics (Figure 11B): Critically, the time delays permit transcription and translation to occur unabatedly until the repression sets in. Specifically, after 3 h, mRNA becomes available, and after an additional 4 h, proteins emerge, including AhRR, which quickly binds to ARNT and begins repressing transcription, resulting soon after in decreased protein production. If the ligand is available beyond time t = 12, the production of protein oscillates, and as soon as the ligand is no longer present, the system returns to the original steady state (not shown). The results for an ODE model without delays and stochasticity are essentially the same as in Figure 11A, and a delay differential equation model, ignoring stochasticity, produces more or less the same results as in Figure 11B. The combination of delays and stochasticity is difficult to capture with differential equations, but it is easily implemented in dBST Figure 11C.

Discussion

Modeling approaches utilizing the framework of Biochemical Systems Theory (BST) have proven powerful in biomedical systems analysis for over seven decades (Savageau, 1969a; Savageau, 1969b; Savageau, 1970). Our goal in the present article was a demonstration that discrete BST (“dBST”) models are noteworthy alternatives to ODE-based BST systems and that they can shed light on complex biomedical phenomena in a similar manner. Discrete dBST models are arguably also more intuitive to newcomers coming from the field of biology, for whom differential equations are often an obscure and dreaded domain of insider mathematics.

The proposal of using dBST is most certainly not a call for abandoning systems of ODEs in biomedical models. ODEs have proven immensely beneficial in all of science, and biomedical applications are no exception. Nonetheless, there are situations that are difficult to align with the concept of instantaneous change. Examples include genuinely discrete events, delays, and stochastic phenomena affecting the phenomenon under study. For instance, we showed elsewhere, in the context of red blood cell death during malarial anemia, that the precise dynamics of blood infections is very difficult to capture with ODEs, but straightforward to implement in a discrete-recursive model (Fonseca and Voit, 2015; Fonseca et al., 2016). Similarly, we demonstrate here and in the Supplementary Data S2 that delays and internal or external stochastic influences affecting a dynamical system are often more easily incorporated into discrete rather than differential equations.

Many of the advantages of BST as a tool for model selection and analysis translate directly into its discrete analog, dBST. Whereas it is generally difficult to choose the most appropriate mathematical formats for representing ill-characterized phenomena a priori, BST and dBST offer guidance at the very beginning of the modeling process, where it is most urgently needed. At the very least, the use of power-law functions offers a viable, unbiased starting point. The power-law format used in BST and dBST is no panacea, but it is a local approximation of mathematically guaranteed quality that typically has a wider range of validity than linear formulations and, embedded into ODEs, is provenly rich enough to permit the inclusion of any differentiable nonlinearities (Voit and Savageau, 1986; Savageau and Voit, 1987).

The use of dBST instead of BST does not create practical design or implementation problems per se, and paradigmatic nonlinearities, such as limit cycles and chaos, can be captured in dBST, as we demonstrated here. If the goal of a dBST model is to mimic a corresponding ODE system as closely as possible, a small step size may have to be chosen. For instance, in the example of limit cycles, a larger step size retained the basic structure and shape of the limit cycle system, but the numerical features were clearly affected. However, the typical task in practical applications is not to create an analog of an ODE system but to convert observed data, together with contextual information, into a computable structure. This inference process is actually simpler in dBST than BST, as most biomedical phenomena are naturally discrete and the determination of optimal parameter values does not require the estimation of slopes.

We demonstrated the ease of designing a dBST model with several small examples and with a moderately complex signal transduction system that triggers changes in gene expression following an exposure to specific toxicants or internal ligands. This phenomenon is difficult to capture with an ODE model because it is critically affected by substantial time delays, which are comingled with the well-known stochastic nature of gene transcription and translation. Our analysis makes it evident that these aspects must not be ignored lest erroneous results are obtained. It also shows how straightforward it is to incorporate these aspects into a dBST model.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/LBSA-VoitLab/dBST.

Author Contributions

Conceived the study: EV; executed all computational analyses: EV and DO; wrote and edited the article: EV and DO.

Funding

This work was supported by the NIEHS: NIH-2P30ES019776-05 (PI: Carmen Marsit); Salary support for EV and DO; Georgia Research Alliance (PI: EV) (general research support; no grant number). The funding agencies are not responsible for the content of this article.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors greatly appreciate constructive feedback from Jacob Davis and Carla Kumbale.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2022.874669/full#supplementary-material

References

Brinkmann, V., Ale-Agha, N., Haendeler, J., and Ventura, N. (2019). The Aryl Hydrocarbon Receptor (AhR) in the Aging Process: Another Puzzling Role for This Highly Conserved Transcription Factor. Front. Physiol. 10, 1561. doi:10.3389/fphys.2019.01561

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlton, J. G., Jones, H., and Eggert, U. S. (2020). Membrane and Organelle Dynamics during Cell Division. Nat. Rev. Mol. Cel Biol. 21, 151–166. doi:10.1038/s41580-019-0208-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Chou, I.-C., and Voit, E. O. (2009). Recent Developments in Parameter Estimation and Structure Identification of Biochemical and Genomic Systems. Math. biosciences 219, 57–83. doi:10.1016/j.mbs.2009.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Claude, D. (1995). Shift of a Limit Cycle in Biology: From Pathological to Physiological Homeostasia*. Chaos 5, 162–166. doi:10.1063/1.166099

PubMed Abstract | CrossRef Full Text | Google Scholar

Cornish-Bowden, A. (2012). Fundamentals of Enzyme Kinetics. Chichester, UK: J. Wiley & Sons, 510.

Google Scholar

Dam, P., Rodriguez-R, L. M., Luo, C., Hatt, J., Tsementzi, D., Konstantinidis, K. T., et al. (2020). Model-based Comparisons of the Abundance Dynamics of Bacterial Communities in Two Lakes. Sci. Rep. 10, 2423. doi:10.1038/s41598-020-58769-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, B. R., Karchner, S. I., Allan, L. L., Pollenz, R. S., Tanguay, R. L., Jenny, M. J., et al. (2008). Repression of Aryl Hydrocarbon Receptor (AHR) Signaling by AHR Repressor: Role of DNA Binding and Competition for AHR Nuclear Translocator. Mol. Pharmacol. 73, 387–398. doi:10.1124/mol.107.040204

PubMed Abstract | CrossRef Full Text | Google Scholar

Fonseca, L. L., Alezi, H. S., Barnwell, J. W., Galinski, M. R., and Voit, E. O. (2016). Quantifying the Removal of Red Blood Cells in Macaca mulatta during a Plasmodium Coatneyi Infection. Malar. J. 15, 410. doi:10.1186/s12936-016-1465-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Fonseca, L. L., and Voit, E. O. (2015). Comparison of Mathematical Frameworks for Modeling Erythropoiesis in the Context of Malaria Infection. Math. biosciences 270, 224–236. doi:10.1016/j.mbs.2015.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Fournier, M. L., Paulson, A., Pavelka, N., Mosley, A. L., Gaudenz, K., Bradford, W. D., et al. (2010). Delayed Correlation of mRNA and Protein Expression in Rapamycin-Treated Cells and a Role for Ggc1 in Cellular Sensitivity to Rapamycin. Mol. Cell Proteomics 9, 271–284. doi:10.1074/mcp.m900415-mcp200

PubMed Abstract | CrossRef Full Text | Google Scholar

Franco, R. S. (2012). Measurement of Red Cell Lifespan and Aging. Transfus. Med. Hemother 39, 302–307. doi:10.1159/000342232

PubMed Abstract | CrossRef Full Text | Google Scholar

Gábor, A., and Banga, J. R. (2015). Robust and Efficient Parameter Estimation in Dynamic Models of Biological Systems. BMC Syst. Biol. 9, 74. doi:10.1186/s12918-015-0219-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Gedeon, T., and Bokes, P. (2012). Delayed Protein Synthesis Reduces the Correlation between mRNA and Protein Fluctuations. Biophysical J. 103, 377–385. doi:10.1016/j.bpj.2012.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Gennemark, P., and Wedelin, D. (2009). Benchmarks for Identification of Ordinary Differential Equations from Time Series Data. Bioinformatics 25, 780–786. doi:10.1093/bioinformatics/btp050

PubMed Abstract | CrossRef Full Text | Google Scholar

Gennemark, P., and Wedelin, D. (2007). Efficient Algorithms for Ordinary Differential Equation Model Identification of Biological Systems. IET Syst. Biol. 1, 120–129. doi:10.1049/iet-syb:20050098

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, A., Galinski, M. R., and Voit, E. O. (2022). Dynamic Control Balancing Cell Proliferation and Inflammation Is Crucial for an Effective Immune Response to Malaria. Front. Mol. Biosci. 8, 800721. doi:10.3389/fmolb.2021.800721

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, G., and Elferink, C. J. (2012). A Novel Nonconsensus Xenobiotic Response Element Capable of Mediating Aryl Hydrocarbon Receptor-dependent Gene Expression. Mol. Pharmacol. 81, 338–347. doi:10.1124/mol.111.075952

PubMed Abstract | CrossRef Full Text | Google Scholar

Julliard, W., Fechner, J. H., and Mezrich, J. D. (2014). The Aryl Hydrocarbon Receptor Meets Immunology: Friend or Foe? A Little of Both. Front. Immunol. 5, 458. doi:10.3389/fimmu.2014.00458

PubMed Abstract | CrossRef Full Text | Google Scholar

Keshet, L. (2005). Mathematical Models in Biology. Philadelphia: Society of Industrial and Applied Mathematics.

Google Scholar

Klipp, E., Liebermeister, W., Wierling, C., Kowald, A., Lehrach, H., and Herwig, R. (2016). Systems Biology: A Textbook. Weinheim, Germany: Wiley VCH Verlag.

Google Scholar

Knowles, I., and Renka, R. J. (2014). Methods for Numerical Differentiation of Noisy Data. Electr. J. Diff. Equations 21, 35–246.

Google Scholar

Koussounadis, A., Langdon, S. P., Um, I. H., Harrison, D. J., and Smith, V. A. (2015). Relationship between Differentially Expressed mRNA and mRNA-protein Correlations in a Xenograft Model System. Sci. Rep. 5, 10775. doi:10.1038/srep10775

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, D. C. (1991). “A Qualitative Analysis of S-Systems: Hopf Bifurcations,” in Canonical Nonlinear Modeling: S-System Approach to Understanding Complexity. Editor E. O. Voit (New York, NY: Van Nostrand Reinhold), 304–344.

Google Scholar

Lewis, M. (2005). The Lac Repressor. Comptes Rendus Biologies 328, 521–548. doi:10.1016/j.crvi.2005.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Beyer, A., and Aebersold, R. (2016). On the Dependency of Cellular Protein Levels on mRNA Abundance. Cell 165, 535–550. doi:10.1016/j.cell.2016.03.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Lorenz, E. N. (1963). Deterministic Nonperiodic Flow. J. Atmos. Sci. 20, 130–141. doi:10.1175/1520-0469(1963)020<0130:dnf>2.0.co;2

CrossRef Full Text | Google Scholar

Lotka, A. J. (1956). Elements of Mathematical Biology. New York: Dover Publications.

Google Scholar

May, R. M. (1973). Stability and Complexity in Model Ecosystems. New Jersey, United States: Princeton University Press.

Google Scholar

McAdams, H. H., and Arkin, A. (1997). Stochastic Mechanisms in Gene Expression. Proc. Natl. Acad. Sci. U.S.A. 94, 814–819. doi:10.1073/pnas.94.3.814

PubMed Abstract | CrossRef Full Text | Google Scholar

Michaelis, L., and Menten, M. L. (1913). Die Kinetik der Invertinwirkung. Biochemische Z. 49, 333–369.

Google Scholar

Nelder, J. A., and Mead, R. (1965). A Simplex Method for Function Minimization. Comput. J. 7, 308–313. doi:10.1093/comjnl/7.4.308

CrossRef Full Text | Google Scholar

R Core Team, R. (2018). A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Google Scholar

Raj, A., and van Oudenaarden, A. (2008). Nature, Nurture, or Chance: Stochastic Gene Expression and its Consequences. Cell 135, 216–226. doi:10.1016/j.cell.2008.09.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Saltzman, B. (1962). Finite Amplitude Free Convection as an Initial Value Problem-I. J. Atmos. Sci. 19, 329–341. doi:10.1175/1520-0469(1962)019<0329:fafcaa>2.0.co;2

CrossRef Full Text | Google Scholar

Savageau, M. A. (1970). Biochemical Systems Analysis. 3. Dynamic Solutions Using a Power-Law Approximation. J. Theor. Biol. 26, 215–226. doi:10.1016/s0022-5193(70)80013-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Savageau, M. A. (1976). Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology. Reading, MA: Addison-Wesley Pub. Co. Advanced Book Program. (reprinted 2009).

Google Scholar

Savageau, M. A. (1969). Biochemical Systems Analysis. II. The Steady-State Solutions for an N-Pool System Using a Power-Law Approximation. J. Theor. Biol. 25, 370–379. doi:10.1016/s0022-5193(69)80027-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Savageau, M. A. (1969). Biochemical Systems Analysis. I. Some Mathematical Properties of the Rate Law for the Component Enzymatic Reactions. J. Theor. Biol. 25, 365–369. doi:10.1016/s0022-5193(69)80026-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Savageau, M. A. (2001). Design Principles for Elementary Gene Circuits: Elements, Methods, and Examples. Chaos 11, 142–159. doi:10.1063/1.1349892

PubMed Abstract | CrossRef Full Text | Google Scholar

Savageau, M. A., and Voit, E. O. (1987). Recasting Nonlinear Differential Equations as S-Systems: A Canonical Nonlinear Form. Math. Biosciences 87, 83–115. doi:10.1016/0025-5564(87)90035-6

CrossRef Full Text | Google Scholar

Schafer, K. A. (1998). The Cell Cycle: a Review. Vet. Pathol. 35, 461–478. doi:10.1177/030098589803500601

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, T. W., Budinsky, R. A., and Rowlands, J. C. (2015). A Model for Aryl Hydrocarbon Receptor-Activated Gene Expression Shows Potency and Efficacy Changes and Predicts Squelching Due to Competition for Transcription Co-activators. PLoS One 10, e0127952. doi:10.1371/journal.pone.0127952

PubMed Abstract | CrossRef Full Text | Google Scholar

Soetaert, K., Petzoldt, T., and Setzer, R. W. (2010). Solving Differential Equations in R: Package deSolve. J. Stat. Softw. 33, 1–25. doi:10.18637/jss.v033.i09

PubMed Abstract | CrossRef Full Text | Google Scholar

Spence, K. F., Decker, S., and Sell, K. W. (1970). Bursting Atlantal Fracture Associated with Rupture of the Transverse Ligament. J. Bone Jt. Surg. 52, 543–549. doi:10.2106/00004623-197052030-00013

CrossRef Full Text | Google Scholar

Stein, R. R., Bucci, V., Toussaint, N. C., Buffie, C. G., Rätsch, G., Pamer, E. G., et al. (2013). Ecological Modeling from Time-Series Inference: Insight into Dynamics and Stability of Intestinal Microbiota. Plos Comput. Biol. 9, e1003388. doi:10.1371/journal.pcbi.1003388

PubMed Abstract | CrossRef Full Text | Google Scholar

Stevens, E. A., Mezrich, J. D., and Bradfield, C. A. (2009). The Aryl Hydrocarbon Receptor: a Perspective on Potential Roles in the Immune System. Immunology 127, 299–311. doi:10.1111/j.1365-2567.2009.03054.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Stockinger, B., Meglio, P. D., Gialitakis, M., and Duarte, J. H. (2014). The Aryl Hydrocarbon Receptor: Multitasking in the Immune System. Annu. Rev. Immunol. 32, 403–432. doi:10.1146/annurev-immunol-032713-120245

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres, N. V., and Voit, E. O. (2002). Pathway Analysis and Optimization in Metabolic Engineering. Cambridge, U.K.: Cambridge University Press.

Google Scholar

Varah, J. M. (1982). A Spline Least Squares Method for Numerical Parameter Estimation in Differential Equations. SIAM J. Sci. Stat. Comput. 3, 28–46. doi:10.1137/0903003

CrossRef Full Text | Google Scholar

Voit, E. O. (2017). A First Course in Systems Biology. 2nd Ed. New York, NY: Garland Science.

Google Scholar

Voit, E. O. (2000). Computational Analysis of Biochemical Systems: A Practical Guide for Biochemists and Molecular Biologists. Cambridge, UK: Cambridge University Press.

Google Scholar

Voit, E. O., Martens, H. A., and Omholt, S. W. (2015). 150 Years of the Mass Action Law. Plos Comput. Biol. 11, e1004012. doi:10.1371/journal.pcbi.1004012

PubMed Abstract | CrossRef Full Text | Google Scholar

Voit, E. O., and Savageau, M. A. (1982). Power-law Approach to Modeling Biological Systems; III. Methods of Analysis. J. Ferment. Technol. 60, 223–241.

Google Scholar

Voit, E. O. (2020). Systems Biology. Oxford, U.K.: Oxford University Press.

Google Scholar

Voit, E. O. (2017). The Best Models of Metabolism. Wiley Interdiscip. Rev. Syst. Biol. Med. 9, e1391. doi:10.1002/wsbm.1391

PubMed Abstract | CrossRef Full Text | Google Scholar

Voit, E. O., and Almeida, J. (2004). Decoupling Dynamical Systems for Pathway Identification from Metabolic Profiles. Bioinformatics 20, 1670–1681. doi:10.1093/bioinformatics/bth140

PubMed Abstract | CrossRef Full Text | Google Scholar

Voit, E. O. (2013). Biochemical Systems Theory: A Review. ISRN Biomathematics 2013, 1–53. doi:10.1155/2013/897658

CrossRef Full Text | Google Scholar

Voit, E. O. (2008). Modelling Metabolic Networks Using Power-Laws and S-Systems. Essays Biochem. 45, 29–40. doi:10.1042/bse0450029

PubMed Abstract | CrossRef Full Text | Google Scholar

Voit, E. O., and Savageau, M. A. (1986). Equivalence between S-Systems and Volterra Systems. Math. Biosciences 78, 47–55. doi:10.1016/0025-5564(86)90030-1

CrossRef Full Text | Google Scholar

Volterra, V. (1926). Variazioni e fluttuazioni del numero d'individui in specie animali conviventi. Mem. R. Accad. Dei Lincei. 2, 31–113.

Google Scholar

Yin, W., and Voit, E. O. (2008). Construction and Customization of Stable Oscillation Models in Biology. J. Biol. Syst. 16, 463–478. doi:10.1142/s0218339008002502

CrossRef Full Text | Google Scholar

Zhu, K., Meng, Q., Zhang, Z., Yi, T., He, Y., Zheng, J., et al. (2019). Aryl Hydrocarbon Receptor Pathway: Role, Regulation and Intervention in Atherosclerosis Therapy (Review). Mol. Med. Rep. 20, 4763–4773. doi:10.3892/mmr.2019.10748

PubMed Abstract | CrossRef Full Text | Google Scholar

Zudaire, E., Cuesta, N., Murty, V., Woodson, K., Adams, L., Gonzalez, N., et al. (2008). The Aryl Hydrocarbon Receptor Repressor Is a Putative Tumor Suppressor Gene in Multiple Human Cancers. J. Clin. Invest. 118, 640–650. doi:10.1172/JCI30024

PubMed Abstract | CrossRef Full Text | Google Scholar