Bulletin of Mathematical Biology

Enciso, German; Kim, Jinsu

doi: 10.1007/s11538-019-00575-3pmid: 30734148

We provide a short review of stochastic modeling in chemical reaction networks for mathematical and quantitative biologists. We use as case studies two publications appearing in this issue of the Bulletin, on the modeling of quasi-steady-state approximations and cell polarity. Reasons for the relevance of stochastic modeling are described along with some common differences between stochastic and deterministic models.

Xu, Bin; Kang, Hye-Won; Jilkine, Alexandra

doi: 10.1007/s11538-019-00573-5pmid: 30756233

Oscillations occur in a wide variety of essential cellular processes, such as cell cycle progression, circadian clocks and calcium signaling in response to stimuli. It remains unclear how intrinsic stochasticity can influence these oscillatory systems. Here, we focus on oscillations of Cdc42 GTPase in fission yeast. We extend our previous deterministic model by Xu and Jilkine to construct a stochastic model, focusing on the fast diffusion case. We use SSA (Gillespie’s algorithm) to numerically explore the low copy number regime in this model, and use analytical techniques to study the long-time behavior of the stochastic model and compare it to the equilibria of its deterministic counterpart. Numerical solutions suggest noisy limit cycles exist in the parameter regime in which the deterministic system converges to a stable limit cycle, and quasi-cycles exist in the parameter regime where the deterministic model has a damped oscillation. Near an infinite period bifurcation point, the deterministic model has a sustained oscillation, while stochastic trajectories start with an oscillatory mode and tend to approach deterministic steady states. In the low copy number regime, metastable transitions from oscillatory to steady behavior occur in the stochastic model. Our work contributes to the understanding of how stochastic chemical kinetics can affect a finite-dimensional dynamical system, and destabilize a deterministic steady state leading to oscillations.

Kang, Hye-Won; KhudaBukhsh, Wasiur; Koeppl, Heinz; Rempała, Grzegorz

doi: 10.1007/s11538-019-00574-4pmid: 30756234

The paper outlines a general approach to deriving quasi-steady-state approximations (QSSAs) of the stochastic reaction networks describing the Michaelis–Menten enzyme kinetics. In particular, it explains how different sets of assumptions about chemical species abundance and reaction rates lead to the standard QSSA, the total QSSA, and the reverse QSSA. These three QSSAs have been widely studied in the literature in deterministic ordinary differential equation settings, and several sets of conditions for their validity have been proposed. With the help of the multiscaling techniques introduced in Ball et al. (Ann Appl Probab 16(4):1925–1961, 2006), Kang and Kurtz (Ann Appl Probab 23(2):529–583, 2013), it is seen that the conditions for deterministic QSSAs largely agree (with some exceptions) with the ones for stochastic QSSAs in the large-volume limits. The paper also illustrates how the stochastic QSSA approach may be extended to more complex stochastic kinetic networks like, for instance, the enzyme–substrate–inhibitor system.

Mai, Ali; Sun, Guowei; Wang, Lin

doi: 10.1007/s11538-018-00568-8pmid: 30627954

In this paper, we study a predator–prey system with random predator dispersal over two habitat patches. We show that in most cases the dispersal delay does not affect the stability and instability of the coexistence equilibrium. However, if the mean time that the predator spent in one patch is much shorter than the timescale of reproduction of the prey and is larger than the double mean time of capture of prey, the dispersal delay can induce stability switches such that an otherwise unstable coexistence equilibrium can be stabilized over a finite number of stability intervals.

Asik, Lale; Peace, Angela

doi: 10.1007/s11538-018-00567-9pmid: 30635835

Phosphorus is an essential element for all life forms, and it is also a limiting nutrient in many aquatic ecosystems. To keep track of the mismatch between the grazer’s phosphorus requirement and producer phosphorus content, stoichiometric models have been developed to explicitly incorporate food quality and food quantity. Most stoichiometric models have suggested that the grazer dynamics heavily depends on the producer phosphorus content when the producer has insufficient nutrient content [low phosphorus (P):carbon (C) ratio]. However, recent laboratory experiments have shown that the grazer dynamics are also affected by excess producer nutrient content (extremely high P:C ratio). This phenomenon is known as the “stoichiometric knife edge.” While the Peace et al. (Bull Math Biol 76(9):2175–2197, 2014) model has captured this phenomenon, it does not explicitly track P loading of the aquatic environment. Here, we extend the Peace et al. (2014) model by mechanistically deriving and tracking P loading in order to investigate the growth response of the grazer to the producer of varying P:C ratios. We analyze the dynamics of the system such as boundedness and positivity of the solutions, existence and stability conditions of boundary equilibria. Bifurcation diagram and simulations show that our model behaves qualitatively similar to the Peace et al. (2014) model. The model shows that the fate of the grazer population can be very sensitive to P loading. Furthermore, the structure of our model can easily be extended to incorporate seasonal P loading.

Marques, Joice; Malchow, Horst; Rodrigues, Luiz; Mistro, Diomar

doi: 10.1007/s11538-018-00539-zpmid: 30644066

The effects of demographic and environmental noise on the vital dynamics and spatial pattern formation are studied for a predator–prey system with strong Allee effect in the prey species. Time and space are taken discrete. It is shown that noise can promote extinction depending on the growth and interaction parameters as well as the noise type and amplitude. The extinction risk increases with the noise amplitude; however, the environmental and demographic noise can have different effects on the risk of extinction. In space, the spatial structures obtained are blurred versions of the deterministic ones in most scenarios. In particular, the complex spatial structures that appear in the parameter domains where the deterministic local dynamics leads to extinction are robust to the density-dependent stochastic fluctuations but are disrupted with environmental noise.

Pages, Nathan; Vera-Sigüenza, Elías; Rugis, John; Kirk, Vivien; Yule, David; Sneyd, James

doi: 10.1007/s11538-018-00563-zpmid: 30644065

We have constructed a spatiotemporal model of $$\hbox {Ca}^{2+}$$ Ca 2 + dynamics in parotid acinar cells, based on new data about the distribution of inositol trisphophate receptors (IPR). The model is solved numerically on a mesh reconstructed from images of a cluster of parotid acinar cells. In contrast to our earlier model (Sneyd et al. in J Theor Biol 419:383–393. https://doi.org/10.1016/j.jtbi.2016.04.030 , 2017b), which cannot generate realistic $$\hbox {Ca}^{2+}$$ Ca 2 + oscillations with the new data on IPR distribution, our new model reproduces the $$\hbox {Ca}^{2+}$$ Ca 2 + dynamics observed in parotid acinar cells. This model is then coupled with a fluid secretion model described in detail in a companion paper: A mathematical model of fluid transport in an accurate reconstruction of a parotid acinar cell (Vera-Sigüenza et al. in Bull Math Biol. https://doi.org/10.1007/s11538-018-0534-z , 2018b). Based on the new measurements of IPR distribution, we show that Class I models (where $$\hbox {Ca}^{2+}$$ Ca 2 + oscillations can occur at constant [ $$\hbox {IP}_3$$ IP 3 ]) can produce $$\hbox {Ca}^{2+}$$ Ca 2 + oscillations in parotid acinar cells, whereas Class II models (where [ $$\hbox {IP}_3$$ IP 3 ] needs to oscillate in order to produce $$\hbox {Ca}^{2+}$$ Ca 2 + oscillations) are unlikely to do so. In addition, we demonstrate that coupling fluid flow secretion with the $$\hbox {Ca}^{2+}$$ Ca 2 + signalling model changes the dynamics of the $$\hbox {Ca}^{2+}$$ Ca 2 + oscillations significantly, which indicates that $$\hbox {Ca}^{2+}$$ Ca 2 + dynamics and fluid flow cannot be accurately modelled independently. Further, we determine that an active propagation mechanism based on calcium-induced calcium release channels is needed to propagate the $$\hbox {Ca}^{2+}$$ Ca 2 + wave from the apical region to the basal region of the acinar cell.

Kitagawa, Kosaku; Kuniya, Toshikazu; Nakaoka, Shinji; Asai, Yusuke; Watashi, Koichi; Iwami, Shingo

doi: 10.1007/s11538-018-00564-ypmid: 30644067

Mathematical modeling has revealed the quantitative dynamics during the process of viral infection and evolved into an important tool in modern virology. Coupled with analyses of clinical and experimental data, the widely used basic model of viral dynamics described by ordinary differential equations (ODEs) has been well parameterized. In recent years, age-structured models, called “multiscale model,” formulated by partial differential equations (PDEs) have also been developed and become useful tools for more detailed data analysis. However, in general, PDE models are considerably more difficult to subject to mathematical and numerical analyses. In our recently reported study, we successfully derived a mathematically identical ODE model from a PDE model, which helps to overcome the limitations of the PDE model with regard to clinical data analysis. Here, we derive the basic reproduction number from the identical ODE model and provide insight into the global stability of all possible steady states of the ODE model.

Mortveit, Henning; Pederson, Ryan

doi: 10.1007/s11538-018-00565-xpmid: 30656504

We present mathematical techniques for exhaustive studies of long-term dynamics of asynchronous biological system models. Specifically, we extend the notion of $$\kappa $$ κ -equivalence developed for graph dynamical systems to support systematic analysis of all possible attractor configurations that can be generated when varying the asynchronous update order (Macauley and Mortveit in Nonlinearity 22(2):421, 2009). We extend earlier work by Veliz-Cuba and Stigler (J Comput Biol 18(6):783–794, 2011), Goles et al. (Bull Math Biol 75(6):939–966, 2013), and others by comparing long-term dynamics up to topological conjugation: rather than comparing the exact states and their transitions on attractors, we only compare the attractor structures. In general, obtaining this information is computationally intractable. Here, we adapt and apply combinatorial theory for dynamical systems from Macauley and Mortveit (Proc Am Math Soc 136(12):4157–4165, 2008. https://doi.org/10.1090/S0002-9939-09-09884-0 ; 2009; Electron J Comb 18:197, 2011a; Discret Contin Dyn Syst 4(6):1533–1541, 2011b. https://doi.org/10.3934/dcdss.2011.4.1533 ; Theor Comput Sci 504:26–37, 2013. https://doi.org/10.1016/j.tcs.2012.09.015 ; in: Isokawa T, Imai K, Matsui N, Peper F, Umeo H (eds) Cellular automata and discrete complex systems, 2014. https://doi.org/10.1007/978-3-319-18812-6_6 ) to develop computational methods that greatly reduce this computational cost. We give a detailed algorithm and apply it to (i) the lac operon model for Escherichia coli proposed by Veliz-Cuba and Stigler (2011), and (ii) the regulatory network involved in the control of the cell cycle and cell differentiation in the Caenorhabditis elegans vulva precursor cells proposed by Weinstein et al. (BMC Bioinform 16(1):1, 2015). In both cases, we uncover all possible limit cycle structures for these networks under sequential updates. Specifically, for the lac operon model, rather than examining all  $$10! > 3.6 \cdot 10^6$$ 10 ! > 3.6 · 10 6 sequential update orders, we demonstrate that it is sufficient to consider 344 representative update orders, and, more notably, that these 344 representatives give rise to 4 distinct attractor structures. A similar analysis performed for the C. elegans model demonstrates that it has precisely 125 distinct attractor structures. We conclude with observations on the variety and distribution of the models’ attractor structures and use the results to discuss their robustness.

Dougoud, Michaël; Mazza, Christian; Schwaller, Beat; Pecze, László

doi: 10.1007/s11538-019-00569-1pmid: 30689102

Here, we present a theoretical investigation with potential insights on developmental mechanisms. Three biological factors, consisting of two diffusing factors and a cell-autonomous immobile transcription factor are combined with different feedback mechanisms. This results in four different situations or fur patterns. Two of them reproduce classical Turing patterns: (1) regularly spaced spots, (2) labyrinth patterns or straight lines with an initial slope in the activation of the transcription factor. The third situation does not lead to patterns, but results in different homogeneous color tones. Finally, the fourth one sheds new light on the possible mechanisms leading to the formation of piebald patterns exemplified by the random patterns on the fur of some cows’ strains and Dalmatian dogs. Piebaldism is usually manifested as white areas of fur, hair, or skin due to the absence of pigment-producing cells in those regions. The distribution of the white and colored zones does not reflect the classical Turing patterns. We demonstrate that these piebald patterns are of transient nature, developing from random initial conditions and relying on a system’s bistability. We show numerically that the presence of a cell-autonomous factor not only expands the range of reaction diffusion parameters in which a pattern may arise, but also extends the pattern-forming abilities of the reaction–diffusion equations.