A model of auto immune response
 James K. Peterson^{1}Email author,
 Alison M. Kesson^{2, 3} and
 Nicholas J. C. King^{4}
DOI: 10.1186/s128650170208x
© The Author(s) 2017
Published: 21 June 2017
Abstract
Background
In this work, we develop a theoretical model of an auto immune response. This is based on modifications of standard second messenger trigger models using both signalling pathways and diffusion and a macro level dynamic systems approximation to the response of a triggering agent such as a virus, bacteria or environmental toxin.
Results
We show that there, in general, will be self damage effects whenever the triggering agent’s effect on the host can be separated into two distinct classes of cell populations.
In each population, the trigger acts differently and this behavior is mediated by the nonlinear interactions between two signalling agents.
Conclusion
If these interactions satisfy certain critical assumptions this will lead to collateral damage. If the initial triggering agent’s action involves any critical host cell population whose loss can lead to serious host health issues, then there is a much increased probability of host death.
Our model also shows that if the nonlinear interaction assumptions are satisfied, there is a reasonable expectation of oscillatory behavior in host health; i.e. periods of remission.
Keywords
Second messenger models Abstract triggering agent Signalling agent mediation Host self damage and death Oscillation in health levels and remission Auto immune responsesBackground
In [1, 2] we explore a micro level simulation model of a single host’s response to varying levels of West Nile Virus (WNV) infection. In that infection, there is a substantial self damage component and in those papers, we show that this is probably due to the way that the virus infects two cell populations differently. This difference, which involves an larger upregulation of MHC1 sites on the surface of nondividing infected cells over dividing infected cells, is critical in establishing a self damage or collateral damage response. In [3], we develop a macro level model of the nonlinear interactions between two critical signalling agents that mediate the interaction between these two sets of cell populations. In the case of WNV infection, the two signals are the MHC1 upregulation level of the cell and the free WNV antigen level. This macro model allowed us to predict a host’s health response to varying levels of initial virus dose. Hence, we could begin to understand the oscillations in collateral damage and host health that lead to the survival data we see in WNV infections. The more general musings of [3] will now be extended to the setting of auto immune interactions in general. The derivations here use standard ideas from advanced calculus and differential equations.
The CMN model
We assume we have a large population of cells T which consists of cells which are infected or altered in two distinct ways by a trigger V. based on signals I, J and K. These two distinct populations of cells will be labeled M and N. There are also non infected cells, H and non infected cells which will be removed due to auto immune action which we call C, for collateral damage. We will be using the same approach to studying nonlinear interactions that was used in [3].
There are then three nonlinear interaction functions F _{1},F _{2} and F _{3} because we know C,M and N depend on each other’s levels in very complicated ways. Usually, we assume the initial trigger dose V _{ 0 } gives rise to some fraction of infected cells and the effect of the trigger will be different in the two cell populations M and N.
Assumption 1
We assume the number of infected cells is p _{0} V _{ 0 } which is split into p _{1} p _{0} V _{ 0 } in population N and p _{2} p _{0} V _{ 0 } in M, where p _{1}+p _{2}=1.
where we now use a standard subscript scheme to indicate the partials. Now let’s add the signals IFN γ (I), J and K to the mix.
The CDN IJK model
where we now use a standard subscript scheme to indicate the partials. If we hold everything constant except i which we increase to i+δ i, what happens? Increasing the IFN γ level should increase collateral damage. Hence, \(H_{1\boldsymbol {i}}^{o} = +\). What about the other two cell populations, M and N?
Assumption 2
We assume this increase in i has no effect. Hence, \(H_{2\boldsymbol {i}}^{o} = H_{3\boldsymbol {i}}^{o} = 0\).
Now hold everything constant except j and increase j to j+δ j. What happens?
Assumption 3
What happens will depend on what the signals J and K are. We can’t argue much yet. However, we suspect the critical assumption to make is the increase in j causes M and N to decrease. Hence, we assume \(H_{2\boldsymbol {j}}^{o} = \) and \(H_{3\boldsymbol {j}}^{o} = \) for each type of infected cell population.
Next hold everything constant except the signal level k and increase k to k+δ k. What happens?
Assumption 4
We assume our choice of signals gives \(H_{2\boldsymbol {k}}^{o} = +\) and \(H_{3\boldsymbol {k}}^{o} = \).
Now we need to estimate i,j and k.
The IJK model
The analysis of the signs of these partials is next. This is similar to what we did for the previous model. If we hold everything constant except i which we increase to i+δ i, what happens? Increasing the IFN γ level should increase IFN γ.
Assumption 5
For the signals j and k, we assume \(G_{1\boldsymbol {i}}^{\boldsymbol {V_{0}}} = +, G_{2\boldsymbol {i}}^{\boldsymbol {V_{0}}} = 0\). and \(G_{3\boldsymbol {i}}^{\boldsymbol {V_{0}}} = 0\) as well.
Now hold everything constant except j and increase j to j+δ j. What happens?
Assumption 6
We assume \(G_{1\boldsymbol {j}}^{\boldsymbol {V_{0}}} = +\), \(G_{2\boldsymbol {j}}^{\boldsymbol {V_{0}}} = +\) and \(G_{3\boldsymbol {j}}^{\boldsymbol {V_{0}}} = +\).
Now hold everything constant except k and increase k to k+δ k. What happens?
Assumption 7
We assume an increase in k will not effect IFN γ levels, i, and hence \(G_{1\boldsymbol {k}}^{\boldsymbol {V_{0}}} = 0\). An increase in k must imply \(G_{3\boldsymbol {k}}^{\boldsymbol {V_{0}}} = +\). We then assume the signals j and k interact so that \(G_{2\boldsymbol {k}}^{\boldsymbol {V_{0}}} = \).
Oscillations in J and K
where A, R, \(G_{1\boldsymbol {i}}^{o}, \beta \) and δ determine a given model.
A health model
Integration details
The JT is a standard integration by parts.
Building the health model
These parameters depend in complex ways on the initial trigger dose V _{ 0 } and it is very difficult to tease out the details.
Collateral damage
Note that there is variation in the collateral damage due to the nonlinear interactions between the J, and K. This is due to the assumptions we have made on the kinds of nonlinear interactions that occur. The remainder of this paper will necessarily discuss why we think these kind of nonlinearities are possible in a variety of auto immune situations.
but then the real part of the eigenvalues would be negative and we would have to model the exponential term as \(e^{r_{2} \boldsymbol {V_{0}} t}\phantom {\dot {i}\!}\). The induced oscillations would then be damped it would not be possible to see oscillatory grown in the collateral damage function.
Note that the minimal values of health do indeed decline for increasing initial trigger dose (although we see up and down variation if we look at all the minimal values versus initial trigger dose as we show in Fig. 1. What we want to focus on here is that the health at many initial trigger doses oscillates over the time of the simulation. For example, the health plot corresponding to V _{ o }=25 has a very low minimum value. These plots are generated by abstractions of health and collateral damage and so it is not clear at all how to relate them to the health of a real host, but it does suggest that the host health can rebound from a low value. Hence, collateral damage can diminish for a given initial trigger which shows a kind of relapse effect.
Also note, the theoretical model we have built so far generates what we call collateral damage and relates it to general health with no discussion of T Cell recognition of infected cells. For us collateral damage is related to IFN γ signalling which is generated by the lysis of cells in the host. In the West Nile Virus infection studied in [1], the splitting into two separate cell populations due to the WNV antigen causes explicitly changes in the avidity computation of the T Cells that recognize MHC1 complexes on infected cells. These changes lead to an enhanced probability that T Cells will target healthy cells because self proteins become more visible to the adaptive immune system. In [3], we derive the heath and collateral damage model we also develop in this paper by focusing very explicitly on the splitting phenomena. The possibility of self damage is now redirected much more strongly to the splitting into two cell populations which allows for a much more general treatment of self damage. However, it is now time to think more deeply about what the signals J and K would be in specific cases of auto immune response to a trigger.
General trigger models
We now want to think carefully about the signals J and K and the trigger V. Hence, we consider the transcriptional control of a regulatory molecule which can be coopted by an external trigger. A good example is the regulatory molecule N F κ B whose normal action is changed by the WNV antigen. It always plays a role in immune system response but the virus alters what it does in many ways. We have discussed this sort of modeling in [4] for the purpose of approximating computation in excitable neurons and in that paper, we follow the spirit of the semiabstract approach in [5]. We are now going to retask it for our purposes of an auto immune discussion, however, this same way of looking at signalling for the context of altering computation in a cognitive model provides another way of looking at the same idea and we think it is always useful to examine a complicated mechanism from multiple points of view.
Consider a trigger T _{0} which activates a cell surface receptor. Inside the cell, there are always protein kinases that can be activated in a variety of ways. Here we denote a protein kinase by the symbol PK. A common mechanism for such an activation is to add to PK another protein subunit to form the complex . This chain of events looks like this: where CSR denotes a cell surface receptor. then acts to phosphorylate another protein. The cell is filled with large amounts of a transcription factor we will denote by T _{1} and an inhibitory protein for T _{1} we label as \(T_{1}^{\sim }\). This symbol, \(T_{1}^{\sim }\), denotes the complement or anti version of T _{1}. In the cell, T _{1} and \(T_{1}^{\sim }\) are generally joined together in a complex denoted by \(T_{1}/T_{1}^{\sim }\). The addition of \(T_{1}^{\sim }\) to T _{1} prevents T _{1} from being able to access the genome in the nucleus to transcribe its target protein. Using methods similar to those discussed in [4], we find the effects of the trigger T _{0} on the change in protein production \(T_{1}, \delta _{T_{1}}\phantom {\dot {i}\!}\), can be modeled by
for β>>1. From this quick analysis, we can clearly see the potentially explosive effect changes in T _{0} can have on .
for parameters μ _{ M }, μ _{ N }, ε _{ M } and ε _{ N }. Hence, we can estimate a strength level for the trigger effect on each population that leads to increases in fragility.
Computational abstractions
We denote the complex formed by the binding of E _{1} and T _{1} by E _{1}/T _{1}. From Fig. 5, we see that the proportion of T _{1} that binds to the genome (DNA) and initiates protein creation P(T _{1}) is thus s r T _{1}.
The protein created, P(T _{1}), could be many things. Here, let us assume that P(T _{1}) is a protein that binds to the promoter for one of the many proteins that effect MHC1 creation, peptide binding etc. For our purposes, we will call it \(\boldsymbol {\mathcal {Q}}\). Thus, our high level model is \(sE_{1}/rT_{1} \: + \: \text {DNA} \rightarrow \boldsymbol {\mathcal {Q}}\). We therefore increase the concentration of MHCI complexes on the surface of the call, thereby making the cell more visible to the adaptive immune system.
We can model the choice process, r T _{1} or (1−r)B _{1}/T _{1} via a simple sigmoid, \(f(x) = 0.5 (1 + \tanh (\frac {xx_{0}}{g}))\) where the transition rate at x _{0} is \(f^{\prime }(x_{0}) = \frac {1}{2g}\). Hence, the “gain” of the transition can be adjusted by changing the value of g. We assume g is positive. This function can be interpreted as switching from of “low” state 0 to a high state 1 at speed \(\frac {1}{2g}\). Now the function h=r f provides an output in (r,∞). If x is larger than the threshold x _{0}, h rapidly transitions to a high state r. On the other hand, if x is below threshold, the output remains near the low state 0.
Our model of the change in \(\boldsymbol {\mathcal {Q}}\) expression is therefore \(\Delta \: \boldsymbol {\mathcal {Q}}^{max} \: = \: h_{\boldsymbol {\mathcal {Q}}}(\left [P(T_{1})\right ])\).
Our model of the change in maximum P(T _{ N }) expression is therefore \(\Delta \: \boldsymbol {P(T_{N})}^{max} \: = \: h_{\boldsymbol {P(T_{N})}}([P(T_{N})])\).
Diffusion trigger models
where D is called the diffusion constant for this substance. For simplicity this is a one dimensional model where the spatial variable x comes from the line segment [0,L]. For example, instead of a three dimensional cell modeled as a sphere, the cell is modeled as a string of finite length. Hence, stuff can only enter the cell from either the right or the left. The condition \(D \frac {\partial u }{ \partial x } \mid _{0,L} = J_{0,L}\) states that there are conditions on the flux of u through the boundary at x=−0 or x=L that must be satisfied. The term J _{0,L } can be thought of as an injection of current into the cell. In this model, we think of the substance u as being in equilibrium in the cell and the current injection J _{0,L } alters that equilibrium and the diffusion model, when solved, tells us what happens to u due to the sudden current injection at the boundary. The critical review on the control of free calcium in cellular processing in [6] notes the concentration of C a ^{++} in the cell is controlled by the reversible binding of calcium ion to the buffer complexes, B _{1},B _{2} and so forth. There in general are quite a few different buffer complexes which all behave differently. These buffer molecules act as calcium ion sensors that, in a sense, decode the information contained in the calcium ion current injection and then pass on a decision to a target. Many of these targets can be proteins transcribed by accessing the genome. Hence, the P(T _{1}) could be a buffer molecule B _{ j } and so the trigger that causes the calcium current injection into the cell could influence the concentration of a buffer B _{ j } and therefore influence how it itself is decoded. In this situation, the boundary condition J _{0,L } plays the role of the entry calcium current. Such a calcium ion input current through the membrane could be due to membrane depolarization causing an influx of calcium ions through the port or via ligand binding to a receptor which in turn indirectly increases free calcium ion in the cytosol. Such mechanisms involve the interplay between the way the endoplasmic reticulum handles release and uptake and also the storage buffers. This boundary current in general therefore determines the u(t,x) solution through the diffusion equation. The exact nature of this solution is determined by the receptor types, buffers and storage sites in the ER. Differences in the period and magnitude of the calcium current u(t,x) resulting from the input J _{0,L } trigger different second messenger pathways. Hence, there are many possible outcomes due to a given input current J _{0,L }.
We will now modify discussions in [7, 8] and [9] that show us how to model calcium ion movement in the cell to develop a model of trigger movement in and out of the cell populations M and N.
We assume the trigger enters the host and can be sequestered in some form in both M and N cell populations. We also assume the trigger has associated with some sort of diffusion process; letting u(t,x) be the concentration of the trigger in the host at time t and spatial position x, we posit u _{ t }=D _{0} u _{ xx }. Also, note we present our arguments as if the host was one dimensional; i.e., the two cell populations lies along a one dimensional axis measured by the location variable x. This is, of course, very simplistic, but we only want to suggest some functional dependencies here, so it will suffice for our purposes. Let’s assume these two populations use or bind the trigger with binding rate constants \(k_{M}^{+}\) and \(k_{N}^{+}\) and disassociation rate constants \(k_{M}^{} \) and \(k_{N}^{}\). Let the total number of cells in the host be P. Then the fraction of cells in the M population is \(\frac {\boldsymbol {M}}{P}\) which we call C _{ M } and the fraction in the N population is \(\frac {\boldsymbol {M}}{P}\) which is denoted by C _{ N }. Hence, the number of cells that have not been exposed to the trigger is P−M(t,x)−N(t,x)=F(t,x).
we obtain
Further approximations
Define the new diffusion constant \(\hat {\mathcal {D}}\) by \(\hat {\mathcal {D}} = \frac {D_{0} + D_{M} \gamma _{M} + D_{N} \gamma _{N}} {\Lambda }\). The free trigger dynamics are thus .
Approximations to trigger modification
where \(e_{P(T_{M})}\) is a scaling factor and \(\delta _{P(T_{M})}\) is the change in P(T _{ M }) expression. Our model of the change in maximum P(T _{ M }) expression is therefore \(\Delta \: \boldsymbol {P(T_{M})}^{max} \: = \: h_{\boldsymbol {P(T_{M})}}([P(T_{M})])\).
Hence, our usual trigger solution is \(u(t,x) = \frac {1}{\sqrt {4 \pi \hat {\mathcal {D}} \: t}} \: e^{\frac {x^{2}}{4 \hat {\mathcal {D}} t}}\) which is altered to \(\hat {u}(t,x) = \frac {1}{\sqrt {4 \pi \tilde {\mathcal {D}} \: t}} \: e^{\frac {x^{2}}{4 \tilde {\mathcal {D}} t}}\).
Results and discussion

The crucial assumption here is that the triggering event has an effect on the host that divides into two parts. For a WNV infection, these two cell populations are the dividing and nondividing infected cells, D and N, respectively. But here, we have posited that these two cell populations are given by M and N instead. Hence, any infectious agents or trigger that gives rise to such a split response engenders a similar collateral damage response which we interpret as an autoimmune reaction. We note this give us a significant insight into general autoimmune responses. Note that Fig. 1 shows there is collateral damage that oscillates due to the trigger and we could interpret a downswing in collateral damage as a relapse event.

We have assumed \(G_{2\boldsymbol {j}}^{\boldsymbol {V_{0}}} = +\), \(G_{3\boldsymbol {j}}^{\boldsymbol {V_{0}}} = +\), \(G_{2\boldsymbol {k}}^{\boldsymbol {V_{0}}} = \), \(G_{3\boldsymbol {k}}^{\boldsymbol {V_{0}}} = +\), which then says the coefficient matrix of the linearized upregulation and free antigen model has the form$$\begin{array}{@{}rcl@{}} \left[\begin{array}{cc} G_{2\boldsymbol{j}}^{\boldsymbol{V_{0}}} & G_{2\boldsymbol{k}}^{\boldsymbol{V_{0}}}\\ G_{3\boldsymbol{j}}^{\boldsymbol{V_{0}}} & G_{3\boldsymbol{k}}^{\boldsymbol{V_{0}}} \end{array}\right] &=& \left[\begin{array}{cc} + & \\ + & + \end{array}\right] \end{array} $$This algebraic sign pattern itself can give rise to complex eigenvalues for the linearized nonlinear interaction model and we have not explored this more general problem. Here, we have posited specific relations that give rise to clearcut oscillations. We have assumed \(G_{2\boldsymbol {j}}^{\boldsymbol {V_{0}}} = G_{3\boldsymbol {k}}^{\boldsymbol {V_{0}}} = \alpha ^{\boldsymbol {V_{0}}}\) and \(G_{3\boldsymbol {j}}^{\boldsymbol {V_{0}}} = G_{2\boldsymbol {k}}^{\boldsymbol {V_{0}}} = \beta ^{\boldsymbol {V_{0}}}\) which gives rise to the characteristic coefficient matrix$$\begin{array}{@{}rcl@{}} \left[\begin{array}{cc} G_{2\boldsymbol{j}}^{\boldsymbol{V_{0}}} & G_{3\boldsymbol{j}}^{\boldsymbol{V_{0}}}\\ G_{3\boldsymbol{j}}^{\boldsymbol{V_{0}}} & G_{2\boldsymbol{J}}^{\boldsymbol{V_{0}}} \end{array}\right] &=& \left[\begin{array}{cc} \alpha^{\boldsymbol{V_{0}}} & \beta^{\boldsymbol{V_{0}}}\\ \beta^{\boldsymbol{V_{0}}} & \alpha^{\boldsymbol{V_{0}}} \end{array}\right] \end{array} $$

We analyze a general trigger in Section “General trigger models” and show that we can understand how the signal generates protein transcription changes via equations such as$$\begin{array}{@{}rcl@{}} \delta_{T_{M}} &=& \mu_{M} \: \left(2 \epsilon_{M}\: + \: \epsilon_{M}^{2}\right) \end{array} $$
which shows how changes in the trigger generate changes in another signal via a cascade of protein transcriptions culminating in a change in fragility for M.

Then, in Section “Computational abstractions”, a computational approach using sigmoid activation function h _{ M } further found that if we associate the change in fragility in M due to T _{ M } to be the alteration of a protein complex P(T _{ M }), then the maximum P(T _{ M }) expression is on the order of \(\Delta \: \boldsymbol {P(T_{M})}^{max} \: = \: h_{\boldsymbol {P(T_{M})}}([P(T_{M})])\) where h is a traditional sigmoid response function which switches the protein from a low to a high level.

We can also model the trigger signal in terms of a diffusion process as we did in Section “Approximations to trigger modification”. The trigger effects M and N is complicated ways and in this section, we study what happens if the pathways the trigger operates on generate a change in M itself. We show this in turn alters the diffusion coefficient that controls the trigger dynamics thereby potentially altering all of the trigger pathways.
from which can conclude we probably have oscillations if the algebraic sign of θ _{3j } is opposite to θ _{2k } and if the algebraic signs of θ _{2j } and θ _{3k } match. We also want θ _{2j } positive so we get undamped oscillations. The more exact equality conditions are probably not actually needed although the analysis was easier when we made them.
How do we use this information? Once we identify signals j and k useful for the dynamical model of interest, we need to experimentally estimate \(\frac {\delta \boldsymbol {j}}{ \boldsymbol {j}}\), \(\frac {\delta \boldsymbol {j} }{\boldsymbol {k}}, \frac {\delta \boldsymbol {k}}{\boldsymbol {j}}\) and \(\frac {\delta \boldsymbol {k}}{\boldsymbol {k}}\). This will give us estimates for the algebraic signs. We believe there will be an autoimmune interaction is the sign patterns we have discussed hold.
Conclusion
In conclusion, we have shown that we can build a reasonable model of how a trigger agent such as a virus, a bacteria or an environmental toxin can infect a host’s cell and cause an autoimmune reaction. Part of our model is a macro one and we believe it provides clues as to how general auto immune reactions behave. We have posited that for an infectious agent or trigger to cause oscillations in health it is required that the trigger causes alterations in two distinct cell populations. Then, if the nonlinear interactions between these two populations satisfies the conditions for damped oscillatory response we have mentioned here, we should see oscillations in the host health and collateral damage. Another part of our model is a detailed micro level one which looks at the triggering pathway and examines how the events in that pathway can contribute to the nonlinear interaction assumptions we make. We also discuss in detail how a sudden trigger increase can be modeled in terms of a diffusion based response and how that also can influence the nonlinear interactions we need for an auto immune response.
We sift through Fig. 7 looking for the following: we want \(A_{pq}^{1} > 0\) for undamped oscillations and we want the oscillation conditions: the sign of \(A_{pq}^{3}\) is the opposite of \(A_{pq}^{2}\) and the signs of \(A_{pq}^{1}\) and \(A_{pq}^{4}\) match. If we can find such a pair (p,q) then we have identified two signals for the two populations of M and N type cells that could imply an auto immune response is possible for this trigger. We think this is quite interesting and could help us decide if a medical event should be classified as an auto immune event.
Methods
We consider this work essentially a theoretical model and we hope that we can generate some insights into the many troubling auto immune problems we face.
Declarations
Funding
Publication of this article was funded by National Health and Medical Research Council Project Grants, Numbers, 512413, 1030897, to NJCK and Australian Research Council Discovery Project Number DP0666152 to AK and NJCK.
Availability of data and materials
Not applicable.
Authors’ contributions
All authors were equally responsible for the ideas and development of the work with JP responsible for the mathematical and computational details. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
About this Supplement
This article has been published as part of BMC Immunology Volume 18 Supplement 1, 2017. Systems Immunology & ImmunoInformatics. The full contents of the supplement are available online https://bmcimmunol.biomedcentral.com/articles/supplements/volume18supplement1.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Peterson J, Kesson AM, King NJC. A Simulation For Flavivirus Infection Decoy Responses. Adv. Microbiol. 2015; 5:123–42. doi:10.4236/aim.201552013.
 Peterson J, Kesson AM, King NJC. Using A Collateral Damage Model To Explain Survival Data In West Nile Virus Infections. Adv Microbiol. 2016; 6(2):251–62. doi:10.4236/aim.2016.64025.View ArticleGoogle Scholar
 Peterson J, Kesson AM, King NJC. A Theoretical Model of the West Nile Virus Survival Data. BMC Immunol. 2016;:1–30. Draft.Google Scholar
 Peterson J. Nodal Computation Approximations in Asynchronous Cognitive Models. Comput Cognit Sci. 2015; 1:4. doi:10.1186/s404690150004y.View ArticleGoogle Scholar
 Gerhart J, Kirschner M. Cells, Embryos and Evolution: Towards a Cellular and Developmental Understanding of Phenotypic Variation and Evolutionary Adaptability. New York: Blackwell Science; 1997.Google Scholar
 Carafoli E, Santella L, Branca D, Brini M. Generation, Control and Processing of Cellular Calcium Signals. Crit Rev Biochem Mol Biol. 2001; 36(2):107–260.View ArticlePubMedGoogle Scholar
 Berridge M. Elementary and Global Aspects of Calcium Signalling. J Physiology. 1997; 499:291–306.View ArticleGoogle Scholar
 Wagner, Keizer. Effects of Rapid Buffers on Ca ^{+2} Diffusion and Ca ^{+2} Oscillations. BioPhys J. 1994; 67:447–56.View ArticlePubMedPubMed CentralGoogle Scholar
 Höffer T, Politi A, Heinrich R. Intracellular Ca ^{+2} Wave Propagation Through Gap  Junctional Ca ^{+2} Diffusion: A Theoretical Study. BioPhys J. 2001; 80(1):75–87.View ArticleGoogle Scholar
 Lancaster M, Renner M, Martin C, Wenzel D, Bicknell L, Hurles M, Homfray T, Penninger J, Jackson A, Knoblich J. Cerebral organoids model human brain development and microcephaly. Nature. 2013; 501:373–9.View ArticlePubMedGoogle Scholar