- Email: [email protected]
Proliferation and Competition in Discrete Biological Systems
Bulletin of Mathematical Biology (2003) 65, 375–396 doi:10.1016/S0092-8240(03)00007-7
Proliferation and Competition in Discrete Biological Systems YORAM LOUZOUN∗ Department of Mathematics, Bar Ilan University, Ramat Gan, 52900, Israel E-mail: [email protected]
SORIN SOLOMON Racah Institute for Physics, Hebrew University, Jerusalem, Israel HENRI ATLAN Human Biology Research Centre, Hadassah Hebrew University Hospital, Jerusalem, Israel IRUN R. COHEN Department of Immunology, The Weizmann Institute of Science, Rehvot, Israel We study the emergence of collective spatio-temporal objects in biological systems by representing individually the elementary interactions between their microscopic components. We use the immune system as a prototype for such interactions. The results of this detailed explicit analysis are compared with the traditional procedure of representing the collective dynamics in terms of densities that obey partial differential equations. The simulations show even for very simple elementary reactions the spontaneous emergence of localized complex structures, from microscopic noise. In turn the effective dynamics of these structures affects the average behaviour of the system in a very decisive way: systems which would according to the differential equations approximation die, display in reality a very lively behaviour. As the optimal modelling method we propose a mixture of microscopic ∗ Author to whom correspondence should be addressed. E-mail: [email protected]
0092-8240/03/030375 + 22 $30.00/0 Society for Mathematical Biology.
c 2003 Published by Elsevier Science Ltd on behalf of
376
Y. Louzoun et al.
simulation systems describing each reaction separately, and continuous methods describing the average behaviour of the agents. c 2003 Published by Elsevier Science Ltd on behalf of Society for Mathematical Biology.
1.
I NTRODUCTION
The immune system is a classical example of biological emergence (Cohen, 2000). It is an archetype for biological systems whose behaviour cannot be directly deduced from their molecular and cellular properties (Kauffman, 1995). As in many other biological systems, a large amount of microscopic data has been gathered. This data by itself is insufficient to explain the macroscopic behaviour of the immune system, and complementary levels of mathematical modelling should accompany it. One of the most widely used mathematical tools in biology is population dynamics described by reaction-diffusion (partial) differential equations (Segel, 1991). We analyse the population dynamics of the immune system and show using a new methodology based on a combination between Monte-Carlo simulation and elements from partial differential equations (PDEs) that traditional methods contain flaws, which make their results unreliable. More precisely, we introduce microscopic simulation (MS) (Stauffer, 1999; Levy et al., 2000; Shnerb et al., 2000a; Louzoun et al., 2001) methods, which on the one hand represent directly the basic elements of the immune system but on the other hand allow the identification of the effective macroscopic objects relevant to the collective dynamics of the system. In the MS the agents taking part in the interactions are placed on a lattice and each reaction is enacted in the computer separately. More precisely we compute at each lattice point the probability of each reaction and perform reactions randomly according to these probabilities. The simulation allows us to analyse which parts of the dynamics are sensitive to the stochastic nature of the interactions and treat them appropriately. Many previous systems tried to deal with the discrete and stochastic nature of biological agents, and specifically in the immune systems (Seiden and Celada, 1992; Meier-Schellersheim and Mack, 1999), but they used very complex agents and interactions between the agents. These systems, although more realistic, fail to discern the basic issues emerging directly from the discrete nature of the agents. In this work, we first show that very simple microscopic interactions lead (when studied with MS) to highly nontrivial effects which one would not have guessed using any of the classical methods. We describe the resulting new phenomena and show how they are relevant for various biological applications. In other words we show that systems that were judged too simple to constitute an adaptive network and too linear to present nonhomogeneous solutions in the partial differential equations context, do in fact
Proliferation and Competition in Discrete Biological Systems
377
emerge well-defined/individualized macroscopic objects, which can move around, act in complex ways, and display adaptive behaviour. We present here MS as a generic technology for representing in a direct intuitive way the microscopic knowledge on the immune system and for extracting nontrivial macroscopic information out of them (including the emergence of adaptive behaviour and self-organization). Our results are in stark contrast with the predictions of the classical modelling methods. In Section 2 we study the macroscopic effect of the proliferation of cells induced by a homogeneous signal distribution. We show that the proliferation rate predicted by the MS is many orders of magnitude larger than the one expected from PDE. In particular one obtains proliferation in conditions that the naive continuum limit predicts decay. In Section 3 we analyse the source of the difference between MS and PDE and explain analytically the errors in the traditional description. In Section 4 we introduce additional reactions, which limit the growth of the proliferating entities. We show that in certain conditions the macroscopic inhomogeneity is even increased by the effect of competition. In conclusion, we are opening the way to a wide front re-evaluation of the dynamics that are supposed to emerge out of the known microscopic elementary interactions of biological systems. We provide the main effects to be re-scrutinized and the particular methods to do so. Paradoxically, the new methods that we propose are closer to the intuition of the practitioners in biology. Indeed the objects which MS manipulates are the very same ones they encounter in their experimental work: cells, molecules and so on (rather than abstract approximations such as differential equations).
2.
DYNAMICS OF D ISCRETE P ROLIFERATING AGENTS
Proliferating entities can be described as agents B that duplicate whenever they meet a homogeneously distributed stimulating agent (a catalyst) A. The B agents also have in general a typical death rate µ. Examples of such agents are: in immunology—immune cells reacting to an antigen (Croft, 1991); in ecology— animals proliferating whenever they find food; in metabolic networks—a protein interaction mediated by a catalyst. This system can be described as a two reactions system (Fig. 1). The first reaction is that whenever a B agent meets a stimulating agent A it multiplies with a rate λ. We will denote this reaction as: B + A → B + B + A; λ. In a lattice formalism, this reaction is interpreted as: ‘whenever a B and an A meet on the same lattice point, a new B is created with a probability of λ’. We will
378
Y. Louzoun et al.
Figure 1. In the right-hand drawing, we describe the reactions taking place in a system of proliferating cells. (1) Whenever a B agent meets an A agent, a new B agent is created. A + B → A + B + B, with a rate of λ. (2) Each B agent has a constant probability µ of dying. B → ·. In the left-hand drawing we present a flow chart of the dynamics. Full line represents activation, while dotted line represents destruction or death.
use the same type of notation in the next sections of the manuscript. The second mechanism is the death of a B agent with a probability µ B → ·; µ. The A and the B agents move around randomly with a diffusion rate D (which can be different for B and for A). One can imagine the following gedanken experiment. We let both the B and the A agents diffuse for a long enough time to assure the homogeneity of their distribution. Then we start the B reactions. The naive lore would expect the growth rate of B to be proportional to number of B and A pairs. Thus the number of new B agents produced in a time interval t at a given location (r) would be: B(r) = A(r)B(r)λt.
(1)
The number of Bs dying is proportional to the total B population: B = −µB(r)t.
(2)
Thus the net local change in B is: B = ( A(r)λ − µ)B(r)t + diffusion
(3)
where diffusion represents the difference between the number of B agents coming from neighbouring points and the number of B agents moving to neighbouring points. This equation becomes a PDE (Boyce and DiPrima, 1986) when we take the limit t → 0 ∂B = (Aλ − µ)B + D∇ 2 B. ∂t
(4)
Proliferation and Competition in Discrete Biological Systems
379
The solution of this PDE with homogeneous initial conditions is Boyce (1986)† : B ∼ e( Aλ−µ)t .
(5)
According to (5), the B population will decay if the average growth rate is lower than the decay rate ((Aλ − µ) > 0) (Malthus, 1798; Verhuulst, 1838). On the other hand, if the average growth rate is higher than the death rate (( Aλ − µ) > 0) the B concentration grows with an exponential growth rate ( Aλ−µ). This exponential growth will end when other mechanisms, discussed in the following sections, become important. In order to test the prediction of the PDE, we simulated the system using discrete agents moving on a lattice. We first let all the agents diffuse freely, until they reach a homogeneous distribution. We then start the interactions described in the previous section. We use a regime where the B birth rate (λ A) is much lower than the death rate µ. If one were to believe the PDE one would conclude that the total B population would decrease to zero. However, the result is quite the opposite. The B population increases exponentially, instead of decreasing to zero (Fig. 2). This simulation was performed using both asynchronous and synchronous simulation algorithms, in order to check for artifacts of the simulation, but the results were similar in both cases. We have shown using a field theory analysis that for a spatially limited system the PDE approximation will fail around a threshold d−2 of DλlA +D B ∼ 1, where l is the grid constant and d is the dimension‡ . However for an unlimited two-dimensional system the PDE approach will eventually always fail (Shnerb et al., 2000b).
3.
S OLUTION OF THE PARADOX
The PDE treatment fails in the case of the previous section, since it ignores the discrete distribution of the A agents, and the correlation between the local B and A concentrations. Indeed, if the signalling agent (A) is a random discrete agent, its density is uniform over large scales but locally its distribution is granular. The discreteness of the A agents induces intrinsic microscopic fluctuations, that are not smeared by diffusion. The diffusion of the A agents changes the configuration of these microscopic fluctuations, but their amplitude does not change. In other words, the A agents are homogeneous either when averaged over a spatial neighbourhood at fixed time, or when averaged over a time interval in a single point. However, locally in time and space there will always be fluctuations in A. † If the initial distribution of B is not homogeneous, the precise solution for the dynamics of B is: 2 ˙ c) are set by the boundary and initial conditions. B(t) = d c φ(c)ec x +( Aλ−µ−c )t . The values of φ(
After a long enough time the homogeneous solution (c = 0) will be much bigger than any other solution and we can ignore the effect of the fluctuations in B. In a paradoxical way we will show that the same type of exponentiation for A fluctuations is the reason for the failure of the PDE treatment. ‡ Hereafter we use a grid constant of 1 to ease the notation.
380
Y. Louzoun et al. Comparison between MS and PDE
1012
1010
B population
108 MS PDE
106
104
102
100
0
5
10
15
20
25
Time
Figure 2. Comparison between the MS and PDE description of the system described in Fig. 1. The dashed line represents the total B agent concentration in the MS and the full line represents the differential equation model. One can observe that in the MS the B population increases indefinitely, while in the PDE the B population decreases to zero. The values used for this simulation are: λ = 0.4, A = 0.5, µ = 1, D A = D B = 0.1, λA − µ = −0.8.
These fluctuations change dramatically the evolution of the average B population. Instead of decreasing rapidly to zero it will increase exponentially to infinity. In the PDE treatment, one assumes that the A fluctuations are smoothed due to the diffusion, and that the time limited fluctuations have no effect on the average behaviour of the system. These two assumptions are wrong. The fluctuations in A cannot be smoothed since they are due to the inherent discreteness of the agents. The effects of the fluctuations cannot be averaged since the A and the B agents are strongly correlated. The local dynamics of the Bs at a point x in space is: ˙ B(x, t) = ( A(x, t)λ − µ)B(x, t) + D B ∇ 2 B.
(6)
In the PDE treatment this dynamics is averaged to: ˙ r = λ Ar Br (x) − µ Br + D B ∇ 2 Br
B = (λ A − µ) Br + D B ∇ 2 Br
(7) (8)
where r represents a region over which the distribution A is smooth enough to be treated by continuous methods. Equation (6) transforms into (7) if one can separate A(x)B(x)r into A(x)r B(x)r , since A(x)r = A. Thus in order for equation (7) to hold A and B have to be independent variables (Papoulis, 1984) at least in the small region r. Equation (7) is the very basis of any PDE approach.
Proliferation and Competition in Discrete Biological Systems
381
Yet, it is wrong in the case presented here, since B(x) is high precisely in the sites where A(x) is high. Thus, A and B are in fact strongly correlated. Moreover, since A is redistributed through diffusion which is a random process, its distribution will converge rapidly to a Poisson distribution. The variance of the A distribution is then as large as its average. These large A fluctuations will induce fluctuation in B which are much larger than the average effect. One could argue that the diffusion of both B and A would smear the local fluctuations of the B population, and weaken the effect of the correlation term. We will show using three limiting cases that this is not the case. We first show that our results are true in the case of low A concentration for a wide range of B and A diffusion rates. We then show the same result through a different mechanism in the case of high A concentration, but low A diffusion. Finally we show that even for a reasonably high A diffusion and concentration and high B concentration our results are true, while PDEs fail. 3.1. Single a approximation. Let us consider first the simplest situation of a single A agent jumping randomly (with a diffusion coefficient D A ) between the locations in an infinite d-dimensional space. We shall denote in the following sections a single A agent as a and a single B agent as b. In between a jumps the B density at the a location grows exponentially (Shnerb et al., 2000a) as B(t) ∼ B(0)e(λ−µ−2d D B )t . The estimation is made by neglecting the flow of B’s from a neighbouring site to the a site. This is justified when the B concentration in the neighbouring sites is much lower than on the a site (it is thus an underestimate of the B growth rate). In the same limit, the ratio between the B density at the a location and the B density on a neighbouring site can be estimated as λ/D B . Consequently, each a jump corresponds to a sudden downwards jump by a factor of λ/D B in the B density. As there are on average 2d D A such jumps per unit time, the net effect of proliferation, diffusion and death, gives the B concentration at the a site as a function of time: B(t) = B(0)e(λ−µ−2d D B −2d D A ln(λ/D B ))t .
(9)
The approximation is in good agreement with the simulation shown in Fig. 3. The slope of the B island maximal density, on a log scale, is indeed ln(λ/D B ), the time dependence of the height of the B island in between a jumps is indeed given approximately by an exponent, (λ − µ − 2d D B )t. Consequently the dashed line (in the inset), which represents equation (3), follows closely the actual growth seen in the simulation (full line). The difference between the theory and simulation is mainly due to the case that two or more a jumps follow each other rapidly, in this case the island’s shape does not stabilize before another a jump is made. These rather rare events, modify somewhat the actual result§ . § The above analysis turns void if D = 0, where the spatial dimensions of the island do not grow B
at all. We do not consider this singular case in this paper. In all other cases, our estimate is an underestimate.
382
Y. Louzoun et al.
Figure 3. Single A approximation. The profile of a B island as a function of time as it follows the random motion of a single A agent. The cross-section of the island is taken through the current location of the A agent. The inset shows the time evolution of the height of the B concentration at the point at which A is currently located (solid line). The B colony grows, although the average growth rate over the entire space is negative (A is extremly low since there is only one A agent in the whole simulation space, thus λA − µ ≈ −µ). The dashed line shows the exponential growth with coefficient λ − µ − 2d D B − 2d D A ln(λD B ). The slope between ‘A jumps’ is λ − µ − 2d D B .
One may ask: what is the situation when B colonies arising around a single a are unstable [i.e., when the exponent in equation (9) is negative]. One possibility is that in such a situation the continuum approximation is valid and the B concentration decays to zero. Another possibility is that, although colonies induced by a single a are unstable, cooperative effects give us back the survival feature. 3.2. Low diffusion approximation. Cooperative effects are clearly observed in the case of a very low A diffusion rate. If the first terms of the exponent are negative (λ − µ − 2d D B < 0), while the last term (−2d D A ln(λ/D B )) is small, a very small nonzero A population would reproduce the exponential growth. For such a low diffusion rate, we can assume that the A configuration is fixed over the time B interactions takes place. For a large enough space, the A will have a Poisson distribution. For a large enough volume, for any average A and for any value m there will be regions in space where Ai >= m (where Ai is the density of A agents at the point i). The m . Thus, if we take probability to have m as at some point is p( Ai = m) = e− A A m! the size of our space to be V e A m! A−m , a large number of points will have at
Proliferation and Competition in Discrete Biological Systems
383
Prediction of MS in the regime of low A diffusion rate
10 12 10 11
B population
10 10 10 9 10 8 10 7 10 6 10 5 10 4
0
10
20
30 Time
40
50
60
Figure 4. The dynamics of the B population in the low A diffusion regime. The parameters used for this simulation are: A = 6, λ = 0.1, µ = 1.0, D B = 0.1, D A = 0.001, λA − µ = −0.4, λ − µ = −0.9. Note that since the A distribution is nearly static the exponent of the B dynamics is constant.
least m as. At every such point we can reproduce the ‘single a approximation’ but with an effective growth rate of λm − µ − 2d D B . If we take m to be big enough then the effective growth rate will be positive and the ‘single a approximation’ will predict an exponential growth. The effective A diffusion rate is m D A , since it is enough for one A to leave the point with the high A density for the effect to fail, but we assume that D A is small enough, and that this effect is negligible. Note that the limiting factor on the value of m is the volume of our space and not only the value of A. Thus even when λ A − µ 0 and λ µ one obtains an exponential growth rate (Fig. 4). 3.3. General theory. We have shown that PDEs fail in two simple to understand cases. The first case represents a low A concentration. In this case the fluctuations are important in comparison with the mean. The second case is a low A diffusion rate. In this case although the A fluctuations are small, they are stable for a long time. We are now in a position to test if we can apply these results to a more general case, where A is not small and D A is the same order of magnitude as the other parameters. We will now show that in a very nonintuitive way the PDE analysis fails even in these cases. This emerging phenomena is completely unpredictable from the microscopic point of view (Fig. 5). In order to understand this phenomena let us decouple the A and B dynamics. If at a given time there is a large A aggregate with m as at a given point, the log of
384
Y. Louzoun et al.
Figure 5. B dynamics within the regime of high A density, mid values of A and B diffusion, λ < µ and λA < µ. The upper plots represents snapshots of the log of the B population at the times pointed by the arrows. One must note the large island that keeps growing through all the snapshots, swallows neighbouring islands and eventually contains a large part of the entire B population, while other islands disappear. The probability that such an island should disappear is inversely proportional to its size. Once a large island is established it has a very low probability of disappearing. However in a limited space after a long enough time even this island will disappear. On the other hand for any given time there is a size of the system for which the total B population keeps growing.
the B population at this point grows following ln(B(x)) = (mλ − µ − 2d D B )t.
(10)
If the local A concentration is not constant, then the B population at this point evolves as: t +t ( A(x, t )λ − µ − 2d D B )dt . (11) ln(B(x)) = t
The log of the ratio between the B population at a small distance r from the point x and the B population at x is: ln(B(x, t)) − ln(B(x + r, t)) =
λr . DB
(12)
Proliferation and Competition in Discrete Biological Systems
385
Let us follow the point of maximal B population in this island. We will denote the original point of maximal B population as x. As long as the integrand on the right-hand side of equation (11) is positive the B population will increase in this site. Once this integrand becomes negative the B population will decrease. However if within the island there are other points where this integrand is positive, the B population can increase in these other points. If the distance of these points from the original B maxima is r, then ln(B) at this point will be at most DλB r lower than in the original position. In order for ln(B) to grow it must increase (on average) in its new position more than DλB r. Condition (12) does not involve the B population. It can be stated entirely in the language of the A dynamics as: t +t λr ( A(x, t )λ − µ − 2d D B )dt > (13) DB t t +t µ + 2d D B r dt > (14) A(x, t ) − ⇒Z= λ DB t t is the interval where the integrand in the left-hand side of equation (14) is positive. Note that the integral itself is defined to be always positive. For the DB , and t is a time interval sake of simplicity we can rescale A to a ∗ = A − µ+2d λ ∗ where a (t) ≥ 0. Equation (14) then becomes: tfinal r a ∗ (t)dt > . (15) DB t0 Let us now analyse the average evolution of an island with a radius of R0 . The average value of r is smaller than R0 . We can thus replace r by R0 in order to obtain a sufficient condition for the island growth. We then use the following algorithm. At every time step we look for the point, within the radius R0 , where Z is maximal. This point can change since the value of a ∗ is changing as the A agents diffuse. Let us name this point the focal point. As long as the value of a ∗ is positive in this point, the B population of the island will grow. Eventually the value of a ∗ at the focal point will become negative. We then reset the value of Z to zero at all the points on the grid, and look again for the point of maximal Z . Note that the point that will become the new focal point can have a negative a ∗ value at the time of the jump, so there can be a small ‘wait period’ between the growth at the previous focal point and the new one. In order for the total B population to grow the average value of Z at the focal points should be bigger then DR0B + wait µλ , where wait is the average wait time. We can thus replace the single A condition equation (9) by a ‘single island condition’: tfinal R0 µ ∗ (16) a (t)dt ≥ + wait . DB λ t0
386
Y. Louzoun et al. 80 70 60
50 40 30 20 10 0
0
0.05
0.1
0.15
0.2
0.25 DA
0.3
0.35
0.4
0.45
0.5
Figure 6. Estimates of Z − wait for a two-dimensional grid with A ∗ = A − 3 and an island diameter of 35. The balance is between the jump needed at the end of the growth period and the growth obtained during the growth period. For any value of D B , one can compute the value of D A for which Z − wait ≥ DR0 . Note that the figure represents B the estimate for a given ‘B island’ diameter, however the estimate obtained from a single diameter is close to the result of the simulation describing the full A and B dynamics.
The new condition is size dependent. In a large enough system, islands of all sizes will temporarily occur. If an island is stable at a given radius it will ‘survive’. In order for the B population to grow, it is thus enough for a single radius to be stable. We check the stability of B islands by testing numerically condition (16) for multiple radiuses¶ . Let us take, for example, a case where λ is half of µ, 2d D B is equal to λ and A is one. In other words the B diffusion rate is as high as the reaction rate, the death rate is much bigger than the birth rate and the A population is not small. According to our scaling relation a ∗ = A − 3. We simulated an island diameter of 35 in two and three dimensions. In two dimensions condition (16) holds for an A diffusion rate of 0.15 and a B diffusion rate of 1.0 (Fig. 6). For three dimensions (not shown) the same condition would hold for A diffusion rates of 0.05. This estimate is obviously an underestimate, since we checked a single diameter, and our algorithm is suboptimal. However simulations with the B agents (Fig. 7) indicates that the computed threshold is reasonable for a system of limited size (100 × 100). ¶ The minimal value of the integral can be computed analytically by ignoring the probability that A
agents shall return to the focal point, and assuming only diffusion away from this point, as we have done in the ‘single A approximation’. However the quality of our estimates depends very much on the return of As to their original position.
Proliferation and Competition in Discrete Biological Systems
387
12
10
10
Total B population
10
8
10
DA=0.1 6
10
4
10
2
DA =0.2
10
0
10
0
50
100
150
Time
Figure 7. Average B cell population on the two sides of the threshold predicted by equation (16). One can see that although the threshold was obtained for a single ‘island’ diameter, it is still a reasonable estimate for systems limited to a size of the order of magnitude of this diameter. When one sets the A diffusion to be lower than the threshold the B population grows indefinitely, while above the threshold the B population decreases to zero.
We have previously shown in the context of a renormalization group analysis that in the regime of high A density and low λ the probability of B survival is much higher in two dimensions than in three dimensions (Shnerb et al., 2000b). This difference is due to the low probability of an a to come back to its original location in three dimensions or more. The current analysis provides an insight on the microscopic details behind our previous results.
4.
M ECHANISMS L IMITING THE P OPULATION G ROWTH
The proliferation of the B agents in the models of Sections 2 and 3 must be limited by some saturation mechanism. We have shown that a high death rate (even a death rate higher than the birth rate) is not enough in order to limit the population proliferation. The required mechanisms can be either an external regulator (Jiang and Chess, 2000), or a limit on the resources available to the proliferating agent (Darwin, 1859). 4.1. Local competition. Many of the mechanisms limiting the B proliferation can be described as competition. Competition for a resource means that whenever two B agents need the resource either to survive or to proliferate only one of them will get it and the other will die (or fail to proliferate) (Fig. 8). In ecology
388
Y. Louzoun et al.
Figure 8. Reaction scheme for a system containing local competition. This system has constant proliferation and death rates as in Fig. 1. However it has an extra mechanism of local competition: Whenever two B agents meet there is a constant probability that one of B. them will be destroyed; B + B
competition can be over food or water, while in the immune system competition can be over access to an antigen. This system is an extension of the logistic equation and can be described using a PDE: ˙ = (λ A(r) − µ)B(r) − γ B(r)2 + D B ∇ 2 B(r). B(r)
(17)
For low B values this system will behave similarly to the system described in the previous sections, where we have already shown the failure of the PDEs. One can analyse this system in two cases. If (λ A − µ) > 0 the total population increases at almost any point in space until limited by the saturation, and then . The discreteness of the A induces stabilizes around a steady state of: B = λ A−µ γ only a low level of fluctuations around this steady state. The more interesting case is when (λ A − µ) < 0. We have shown in the previous sections that, even in this case, the point B = 0 is not a stable steady state (in contradiction with the PDE result). The steady state occurs when the B population in the regions of high B population is saturated by the competition term. The only difference between this case and the one presented in the previous sections is the limited size of the high B population zones (Fig. 9). Paradoxically the conditions where PDE and the MS show similar features are seldom realized in nature. They represent extreme cases where a certain agent is able to fill all the available space. In the animal world only a few species enjoy such a status (for example, humans and their parasites). In immunology this regime corresponds to a cell type filling all the blood and immune systems (systemic diseases). Except for these extreme cases the situation where (λ A − µ) < 0 seems to be more appropriate for realistic systems.
Proliferation and Competition in Discrete Biological Systems
389
Figure 9. Comparison between PDE and MS for the system defined by Fig. 8. This system has an average negative growth rate yielding an exponential decay in the PDE description as in the previous section. In the MS we observe a limited growth of the B islands, which is due to the competition between Bs. Thus even regions where the local B concentration would grow are limited to a finite concentration.
4.2. Global competition. The proliferation of the B agent, and the competition between different B agents can occur over very different spatial scales. The proliferation is a local mechanism; the factor leading to competition can be global. In ecological systems the global character of the limiting factor may be related with the fact that the predators move much faster than their prey and that their population is proportional to the total prey population in their hunting region. Therefore, their population constitutes an indirect interaction between distant prey locations. We will show that in such cases the inhomogeneous effects due to discreteness are even more important than in the previous sections. As an aside let us note that this kind of reactions although formally local (a lion can eat a prey only if they meet) is very effective in inducing top down effects; that is the regulation of the local prey population is enforced by their total population. MS can describe in a natural way only local interactions. It is nontrivial to describe with MS factors acting on widely different scales. The solution in this case would be to create a new type of Hybrid Model, which will allow in particular to express the interaction between agents and their averages over various ranges.
390
Y. Louzoun et al.
Figure 10. Global competition reaction scheme. We added to the reaction described in Fig. 1 a global competition interaction. The death probability of each B agent is proportional to the total B agent population.
A simple example for such a system is obtained by modifying the regulation mechanism of the system presented in the previous section (Fig. 10). There the two mechanisms proliferation and competition were local interactions. We can now replace the local competition by a global competition. Even though from the microscopic point of view the nonlocality might look unnatural this is the only way that nonlinearities can arise at the macroscopic level. Indeed an ODE term of the form B2 can appear only as the average over terms of the type B(r) ∗ B. Therefore the modification of equation (6) would be ˙ = λ A(r)B(r) − µB(r) − γ B(r) B + D∇ 2 B(r). B(r)
(18)
That is we replaced the local competition in the system described in Section 4.1 by a global competition term. We find that passing from local to global reactions has crucial implications: by comparing the behaviour of the system equation (18) with the system equation (17) one finds that they are dramatically different. Even in the region where the system equation (17) is described correctly by a PDE, a PDE is still incapable of describing properly the system equation (18). Indeed in system equation (18) (for most parameter sets) the entire space is empty of Bs except for a single island centreed around the point where the A concentration is maximal (Fig. 11). We can understand this system by ignoring (temporarily) the diffusion term and simplifying the equation to : B˙ = A(r)B(r) − B(r) B = B(r)(A(r) − B).
(19)
We rescaled the equation to λ = 1 and included the death term µ into the A density. This can be done with no loss of generality.
Proliferation and Competition in Discrete Biological Systems
391
Figure 11. Comparison between PDE and MS for a system with an average negative growth rate and a global competition term (as defined in Fig. 10). In this model a single large island with a very high B concentration appears. The competition with this island destroys all the other small B islands. The (a) and (b) drawings show the log of the B agents concentration at different times. The (c) drawing shows the comparison between the PDE and the MS.
One can see that in every point that A(r) is lower than Amax the B(r) density will vanish, while around the point rmax defined by A(r) = Amax a region of high B density will appear. Suppose that one starts with all the B(r)s very small; as long as the average concentration of the Bs ( B) is lower than Amax , the local concentration B(r) at the point rmax will rise [according to equation (19)]. This in turn will raise the value of B. Eventually B will reach the value of Amax . What will be the equilibrium value of the B(r)s in every other point? In every point where A(r) < Amax , one will have B(r)( A(r) − B) = B(r)( A(r) − Amax ) < 0, which according to (19) means that the concentration at those points will decay to 0. It is clear therefore that the entire B population responsible for B = Amax is due to a very dense island situated around rmax . The equality B = Amax is not affected by the B diffusion, which only has the effect of spreading the distribution of the B(r)’s in a limited region around rmax . One sees that the dynamics is dominated in this case by the most extreme stochastic local fluctuation in the A(r) distribution, and not by the global or average properties of the A(r)s.
392
Y. Louzoun et al.
Figure 12. The model defined in Fig. 9 in a regime of high A diffusion rate. When the A diffusion rate is very high, a few large B islands are created. As oppossed to Fig. 11, the largest B island does not always dominate, and smaller island temporarily appear and disappear. The various islands emerge an adaptive behaviour: they look for the maxima of the A concentration, split, merge and die. The top windows shows different snapshots of the B distribution. The bottom window is the long-term evolution of this system showing intermittent fluctuations.
In the presence of A diffusion two regimes appear. In most of phase space all the B population is concentrated in a single island, but in a small part of phase space, multiple islands can exist for a long time (Fig. 12). The phase transition occurs when the destruction time of small B islands is longer than the time over which the λA −µ−2d D B , largest island is stable. The total B density is determined by Btot = big isl γ where Abig isl is the A density at the centre of the largest island. The disappearing 1 1 = λ( Abig isl −A . time of smaller island is proportional to λ Asmall isl −µ−2d D B −γ Btot small isl ) The difference between the A density in the largest islands and the next smaller islands is of the order of unities. The largest island will thus lose its advantage over other islands at a rate proportional to D A . Multiple islands can thus exist if DA ∼ DλA > 1, on the other hand if D A is much bigger than λ then λ( Abig isl −Asmall isl ) no island will survive and the results will converge to the PDE estimate. We are thus left with a small part of shape space close to the phase transition in which multiple islands can survive. Within this regime, the diffusion of the A does not
Proliferation and Competition in Discrete Biological Systems
393
lead to smoothing of the B fluctuations. It only leads to a movement of the high B population zone. This zone does not disappear, since its lifespan is much larger than the lifespan of a single B (it is the time it takes to destroy a macroscopically large number of Bs). Within this lifespan a new large fluctuation of A will appear within this zone, and this zone will start to grow again. The b islands population is then greatly affected by the local changes in the a distribution and the entire system undergoes intermittent macroscopic fluctuations [in certain systems such fluctuations can have catastrophic consequences (Louzoun et al., 2002)]. 4.3. Emergence of complex behaviour. In the first sections it has been shown that the discrete character of the elementary agents A leads to unexpectedly rich and complex behaviour of the B population. The complex spatio-temporal emerging structure of the B population cannot be predicted by a PDE approach. In the present chapter we introduced a few types of competition terms. Local competition limits the size and maximum concentration of each macroscopic B island. Rather than interfering with the emergence of macroscopic features, the competition terms turned out to enhance the local and complex character of the B islands. Global competition can endow these macroscopic islands with an emerging social (more precisely anti-social) behaviour. The evolving macroscopic objects now have apparent ‘goals’: seeking food (As), multiplying, dying or destroying the competing islands. The lifespan of the emerging islands is many orders of magnitude longer than the lifespan of their components (which of course remain nonadaptive and totally mechanical during the entire history of the system).
5.
D ISCUSSION
Dynamic biological systems were traditionally described using PDEs. PDEs were invented and proved useful in simple physical systems that fulfilled certain restrictions: • A very large density of interacting particles. • A diffusion rate high enough to smooth local fluctuations. • An immediate result of each interaction, i.e., no delay between the interaction and its result. These assumptions are often false in cellular systems. • The density of interacting particles is often small; although the total number of interacting agents may be very large, the local number of interacting cells can be very low. In order to use a PDE approximation the Poissonian fluctuations in the number of interacting agents should be much smaller than the average number of interacting agents. Especially in the case of autocatalytic systems the density should be high enough that the probability of a large
394
Y. Louzoun et al.
fluctuation should be low. Note that in contradiction with the naive intuition, the probability of larger fluctuations increases with the size of the modelled space. • The diffusion rate can in many cases be low. The diffusion rate should be defined as the ratio between the rate of the random motion of the agents and their interaction radius. Moreover, compartmentalization such as physical objects in ecology or endoplasmic reticulum in intracellular interactions can further limit the diffusion. To summarize, PDEs can be used only if the agents move freely faster than the interaction radius divided by the reaction rate. • Most biological processes require the production of new molecules or the division of cells. Each of these processes can take between hours and days. Thus, there is a delay between an action and its result. This aspect is the most problematic since even minor delays can change the dynamics of the system. A simple example is the loss of stability of the fix point in predator– prey in the presence of a small delay. At the current point we cannot provide a ‘thumb rule’ for when delays can be ignored. We studied in this paper one of the main sources of complexity: the combination of microscopic inhomogeneity and autocatalysis (Shnerb et al., 2000a). In biological systems, there is an inherent source of microscopic inhomogeneity—the discreteness of the entities composing biological systems. Moreover, most biological systems involve proliferation, which is the basic autocatalytic mechanism. We showed that in some very simple cases containing these two elements PDEs fail to describe appropriately the dynamics. This is as well: systems, which were too simple to present complex emergent features in the framework of PDEs turned out to present when studied appropriately a complex behaviour. We used in addition to the straightforward MS intermediate techniques that allow to take explicitly into account the intermediate relevant dynamical scales between the MS and PDEs. In the systems simulated in our paper, autocatalysis results in the flow of information from small-scale to large-scale features; from the microscopic details of the system to its macroscopic behaviour. The rapid evolution of small-scale features to a macroscopic scale affects the distribution of the macroscopic variables. In particular the macroscopic variables which are dominated by the sum of a small number of large values (instead of the sum of a large number of small values) have a very specific (fractal-power law) space–time behaviour. Such effects are enhanced when there is a direct effect of the macroscopic values on the microscopic mechanisms. Many biological systems evolve at time scales much longer than their fundamental microscopic time scales. In PDEs the emergence of long time scales cannot occur naturally, while in MS they emerge spontaneously. In MS, the agents selforganize (Atlan, 1987) in collective structures with emergent adaptive properties and with dynamical space–time scales of their own. The effective dynamics of these emerging objects is much slower than the one of their components [and this fact constitutes their very definition (Solomon, 1995)]. Analysing the system at the
Proliferation and Competition in Discrete Biological Systems
395
level of the emerging objects automatically requires the transition to longer time scales and larger spatial scale. Each level corresponds to a new space–time scale in a hierarchy of complementary representations of the same complex macroscopic system. To summarize, we found in this paper that one of the core mechanisms for the emergence of collective macroscopic objects is the interplay between autocatalysis and fluctuations inherent to the discrete structure of the microscopic components of the system. PDEs or even stochastic PDEs cannot describe this mechanism, since both the emergence and the disappearance of these collective objects are a result of their discrete nature.
R EFERENCES Atlan, H. (1987). Self-creation of meaning. Phys. Scr. 36, 563–580. Boyce, W. E. and R. C. DiPrima (1986). Elementary Differential Equations and Boundary Value Problems, 4th edn, New York: Wiley. Cohen, I. R. (2000). Tending Adam’s Garden: Evolving the Cognitive Immune Self, San Diego, CA: Academic Press. Croft, M. (1991). Activation of naive, memory and effector T cells. Curr. Opin. Immunol. 6, 431–437. Darwin, C. (1859). The Origin of Species, London, U.K.: John Murray. Jiang, H. and L. Chess (2000). The specific regulation of immune responses by CD8+ T cells restricted by the MHC class Ib molecule, Qa-1. Annu. Rev. Immunol. 18, 185–216. Kauffman, A. S. (1995). At Home in the Universe: the Search of the Laws of Complexity, U.K.: Oxford University Press. Levy, M., H. Levy and S. Solomon (2000). Microscopic Simulation of Financial Markets; From Investor Behavior to Market Phenomena, New York: Academic Press. Louzoun, Y., D. Mazurski, J. Goldberg and S. Solomon (2002). The risk of being unfair, the unstable fate of globalisation (submitted). Louzoun, Y., S. Solomon, H. Atlan and I. R. Cohen (2001). The emergence of spatial complexity in the immune system. Physica A 297, 242–252. Malthus, T. R. (1798). An Essay on the Principle of Population printed for J Johnsonin, St Paul’s Churchyard, London (reprinted by Macmillan et Co 1894). Meier-Schellersheim, M. and G. Mack (1999). SIMMUNE, a tool for simulating and analyzing immune system behavior DESY-99-034. Papoulis, A. (1984). Probability, Random Variables, and Stochastic Processes, 2nd edn, New York: McGraw-Hill, pp. 144–145. Segel A. L. (1991). Biological Kinetics, Cambridge Studies in Mathematical Biology 12, Cambridge, U.K. Seiden, P. E. and F. Celada (1992). A model for simulating cognate recognition and response in the immune system. J. Theor. Biol. 158, 329. Shnerb, N., Y. Louzoun, E. Bettelheim and S. Solomon (2000a). The importance of being discrete—life always wins on the surface. Proc. Natl Acad. Sci. 97, 10322. Shnerb, N., E. Bettelheim, Y. Louzoun, O. Agam and S. Solomon (2000b). Adaptation of autocatalytic fluctuations to diffusive noise. Phys. Rev. E 63, 021103.
396
Y. Louzoun et al.
Solomon, S. (1995). The microscopic representation of complex macroscopic phenomena: critical slowing down—a blessing in disguise, in Annual Reviews of Computational Physics II, D. Stauffer (Ed.), World Scientific, pp. 466. Stauffer, D. (1999). Anual Reviews of Computational Physics VII, Singapore: World Scientific Publishing Company. Verhuulst P. F. (1838). Notice sur la Loi que la Population Suit dans son Accroissement; Correspondence Mathematique et Physique, A. Quetelet (Ed.) pp. 113–121; (English translation in Smith, D. and N. Keyfiz (1977). Mathematical Demography, New York: Springer, pp. 333–337). Received 7 May 2002 and accepted 27 November 2002
Proliferation and Competition in Discrete Biological Systems YORAM LOUZOUN∗ Department of Mathematics, Bar Ilan University, Ramat Gan, 52900, Israel E-mail: [email protected]
SORIN SOLOMON Racah Institute for Physics, Hebrew University, Jerusalem, Israel HENRI ATLAN Human Biology Research Centre, Hadassah Hebrew University Hospital, Jerusalem, Israel IRUN R. COHEN Department of Immunology, The Weizmann Institute of Science, Rehvot, Israel We study the emergence of collective spatio-temporal objects in biological systems by representing individually the elementary interactions between their microscopic components. We use the immune system as a prototype for such interactions. The results of this detailed explicit analysis are compared with the traditional procedure of representing the collective dynamics in terms of densities that obey partial differential equations. The simulations show even for very simple elementary reactions the spontaneous emergence of localized complex structures, from microscopic noise. In turn the effective dynamics of these structures affects the average behaviour of the system in a very decisive way: systems which would according to the differential equations approximation die, display in reality a very lively behaviour. As the optimal modelling method we propose a mixture of microscopic ∗ Author to whom correspondence should be addressed. E-mail: [email protected]
0092-8240/03/030375 + 22 $30.00/0 Society for Mathematical Biology.
c 2003 Published by Elsevier Science Ltd on behalf of
376
Y. Louzoun et al.
simulation systems describing each reaction separately, and continuous methods describing the average behaviour of the agents. c 2003 Published by Elsevier Science Ltd on behalf of Society for Mathematical Biology.
1.
I NTRODUCTION
The immune system is a classical example of biological emergence (Cohen, 2000). It is an archetype for biological systems whose behaviour cannot be directly deduced from their molecular and cellular properties (Kauffman, 1995). As in many other biological systems, a large amount of microscopic data has been gathered. This data by itself is insufficient to explain the macroscopic behaviour of the immune system, and complementary levels of mathematical modelling should accompany it. One of the most widely used mathematical tools in biology is population dynamics described by reaction-diffusion (partial) differential equations (Segel, 1991). We analyse the population dynamics of the immune system and show using a new methodology based on a combination between Monte-Carlo simulation and elements from partial differential equations (PDEs) that traditional methods contain flaws, which make their results unreliable. More precisely, we introduce microscopic simulation (MS) (Stauffer, 1999; Levy et al., 2000; Shnerb et al., 2000a; Louzoun et al., 2001) methods, which on the one hand represent directly the basic elements of the immune system but on the other hand allow the identification of the effective macroscopic objects relevant to the collective dynamics of the system. In the MS the agents taking part in the interactions are placed on a lattice and each reaction is enacted in the computer separately. More precisely we compute at each lattice point the probability of each reaction and perform reactions randomly according to these probabilities. The simulation allows us to analyse which parts of the dynamics are sensitive to the stochastic nature of the interactions and treat them appropriately. Many previous systems tried to deal with the discrete and stochastic nature of biological agents, and specifically in the immune systems (Seiden and Celada, 1992; Meier-Schellersheim and Mack, 1999), but they used very complex agents and interactions between the agents. These systems, although more realistic, fail to discern the basic issues emerging directly from the discrete nature of the agents. In this work, we first show that very simple microscopic interactions lead (when studied with MS) to highly nontrivial effects which one would not have guessed using any of the classical methods. We describe the resulting new phenomena and show how they are relevant for various biological applications. In other words we show that systems that were judged too simple to constitute an adaptive network and too linear to present nonhomogeneous solutions in the partial differential equations context, do in fact
Proliferation and Competition in Discrete Biological Systems
377
emerge well-defined/individualized macroscopic objects, which can move around, act in complex ways, and display adaptive behaviour. We present here MS as a generic technology for representing in a direct intuitive way the microscopic knowledge on the immune system and for extracting nontrivial macroscopic information out of them (including the emergence of adaptive behaviour and self-organization). Our results are in stark contrast with the predictions of the classical modelling methods. In Section 2 we study the macroscopic effect of the proliferation of cells induced by a homogeneous signal distribution. We show that the proliferation rate predicted by the MS is many orders of magnitude larger than the one expected from PDE. In particular one obtains proliferation in conditions that the naive continuum limit predicts decay. In Section 3 we analyse the source of the difference between MS and PDE and explain analytically the errors in the traditional description. In Section 4 we introduce additional reactions, which limit the growth of the proliferating entities. We show that in certain conditions the macroscopic inhomogeneity is even increased by the effect of competition. In conclusion, we are opening the way to a wide front re-evaluation of the dynamics that are supposed to emerge out of the known microscopic elementary interactions of biological systems. We provide the main effects to be re-scrutinized and the particular methods to do so. Paradoxically, the new methods that we propose are closer to the intuition of the practitioners in biology. Indeed the objects which MS manipulates are the very same ones they encounter in their experimental work: cells, molecules and so on (rather than abstract approximations such as differential equations).
2.
DYNAMICS OF D ISCRETE P ROLIFERATING AGENTS
Proliferating entities can be described as agents B that duplicate whenever they meet a homogeneously distributed stimulating agent (a catalyst) A. The B agents also have in general a typical death rate µ. Examples of such agents are: in immunology—immune cells reacting to an antigen (Croft, 1991); in ecology— animals proliferating whenever they find food; in metabolic networks—a protein interaction mediated by a catalyst. This system can be described as a two reactions system (Fig. 1). The first reaction is that whenever a B agent meets a stimulating agent A it multiplies with a rate λ. We will denote this reaction as: B + A → B + B + A; λ. In a lattice formalism, this reaction is interpreted as: ‘whenever a B and an A meet on the same lattice point, a new B is created with a probability of λ’. We will
378
Y. Louzoun et al.
Figure 1. In the right-hand drawing, we describe the reactions taking place in a system of proliferating cells. (1) Whenever a B agent meets an A agent, a new B agent is created. A + B → A + B + B, with a rate of λ. (2) Each B agent has a constant probability µ of dying. B → ·. In the left-hand drawing we present a flow chart of the dynamics. Full line represents activation, while dotted line represents destruction or death.
use the same type of notation in the next sections of the manuscript. The second mechanism is the death of a B agent with a probability µ B → ·; µ. The A and the B agents move around randomly with a diffusion rate D (which can be different for B and for A). One can imagine the following gedanken experiment. We let both the B and the A agents diffuse for a long enough time to assure the homogeneity of their distribution. Then we start the B reactions. The naive lore would expect the growth rate of B to be proportional to number of B and A pairs. Thus the number of new B agents produced in a time interval t at a given location (r) would be: B(r) = A(r)B(r)λt.
(1)
The number of Bs dying is proportional to the total B population: B = −µB(r)t.
(2)
Thus the net local change in B is: B = ( A(r)λ − µ)B(r)t + diffusion
(3)
where diffusion represents the difference between the number of B agents coming from neighbouring points and the number of B agents moving to neighbouring points. This equation becomes a PDE (Boyce and DiPrima, 1986) when we take the limit t → 0 ∂B = (Aλ − µ)B + D∇ 2 B. ∂t
(4)
Proliferation and Competition in Discrete Biological Systems
379
The solution of this PDE with homogeneous initial conditions is Boyce (1986)† : B ∼ e( Aλ−µ)t .
(5)
According to (5), the B population will decay if the average growth rate is lower than the decay rate ((Aλ − µ) > 0) (Malthus, 1798; Verhuulst, 1838). On the other hand, if the average growth rate is higher than the death rate (( Aλ − µ) > 0) the B concentration grows with an exponential growth rate ( Aλ−µ). This exponential growth will end when other mechanisms, discussed in the following sections, become important. In order to test the prediction of the PDE, we simulated the system using discrete agents moving on a lattice. We first let all the agents diffuse freely, until they reach a homogeneous distribution. We then start the interactions described in the previous section. We use a regime where the B birth rate (λ A) is much lower than the death rate µ. If one were to believe the PDE one would conclude that the total B population would decrease to zero. However, the result is quite the opposite. The B population increases exponentially, instead of decreasing to zero (Fig. 2). This simulation was performed using both asynchronous and synchronous simulation algorithms, in order to check for artifacts of the simulation, but the results were similar in both cases. We have shown using a field theory analysis that for a spatially limited system the PDE approximation will fail around a threshold d−2 of DλlA +D B ∼ 1, where l is the grid constant and d is the dimension‡ . However for an unlimited two-dimensional system the PDE approach will eventually always fail (Shnerb et al., 2000b).
3.
S OLUTION OF THE PARADOX
The PDE treatment fails in the case of the previous section, since it ignores the discrete distribution of the A agents, and the correlation between the local B and A concentrations. Indeed, if the signalling agent (A) is a random discrete agent, its density is uniform over large scales but locally its distribution is granular. The discreteness of the A agents induces intrinsic microscopic fluctuations, that are not smeared by diffusion. The diffusion of the A agents changes the configuration of these microscopic fluctuations, but their amplitude does not change. In other words, the A agents are homogeneous either when averaged over a spatial neighbourhood at fixed time, or when averaged over a time interval in a single point. However, locally in time and space there will always be fluctuations in A. † If the initial distribution of B is not homogeneous, the precise solution for the dynamics of B is: 2 ˙ c) are set by the boundary and initial conditions. B(t) = d c φ(c)ec x +( Aλ−µ−c )t . The values of φ(
After a long enough time the homogeneous solution (c = 0) will be much bigger than any other solution and we can ignore the effect of the fluctuations in B. In a paradoxical way we will show that the same type of exponentiation for A fluctuations is the reason for the failure of the PDE treatment. ‡ Hereafter we use a grid constant of 1 to ease the notation.
380
Y. Louzoun et al. Comparison between MS and PDE
1012
1010
B population
108 MS PDE
106
104
102
100
0
5
10
15
20
25
Time
Figure 2. Comparison between the MS and PDE description of the system described in Fig. 1. The dashed line represents the total B agent concentration in the MS and the full line represents the differential equation model. One can observe that in the MS the B population increases indefinitely, while in the PDE the B population decreases to zero. The values used for this simulation are: λ = 0.4, A = 0.5, µ = 1, D A = D B = 0.1, λA − µ = −0.8.
These fluctuations change dramatically the evolution of the average B population. Instead of decreasing rapidly to zero it will increase exponentially to infinity. In the PDE treatment, one assumes that the A fluctuations are smoothed due to the diffusion, and that the time limited fluctuations have no effect on the average behaviour of the system. These two assumptions are wrong. The fluctuations in A cannot be smoothed since they are due to the inherent discreteness of the agents. The effects of the fluctuations cannot be averaged since the A and the B agents are strongly correlated. The local dynamics of the Bs at a point x in space is: ˙ B(x, t) = ( A(x, t)λ − µ)B(x, t) + D B ∇ 2 B.
(6)
In the PDE treatment this dynamics is averaged to: ˙ r = λ Ar Br (x) − µ Br + D B ∇ 2 Br
B = (λ A − µ) Br + D B ∇ 2 Br
(7) (8)
where r represents a region over which the distribution A is smooth enough to be treated by continuous methods. Equation (6) transforms into (7) if one can separate A(x)B(x)r into A(x)r B(x)r , since A(x)r = A. Thus in order for equation (7) to hold A and B have to be independent variables (Papoulis, 1984) at least in the small region r. Equation (7) is the very basis of any PDE approach.
Proliferation and Competition in Discrete Biological Systems
381
Yet, it is wrong in the case presented here, since B(x) is high precisely in the sites where A(x) is high. Thus, A and B are in fact strongly correlated. Moreover, since A is redistributed through diffusion which is a random process, its distribution will converge rapidly to a Poisson distribution. The variance of the A distribution is then as large as its average. These large A fluctuations will induce fluctuation in B which are much larger than the average effect. One could argue that the diffusion of both B and A would smear the local fluctuations of the B population, and weaken the effect of the correlation term. We will show using three limiting cases that this is not the case. We first show that our results are true in the case of low A concentration for a wide range of B and A diffusion rates. We then show the same result through a different mechanism in the case of high A concentration, but low A diffusion. Finally we show that even for a reasonably high A diffusion and concentration and high B concentration our results are true, while PDEs fail. 3.1. Single a approximation. Let us consider first the simplest situation of a single A agent jumping randomly (with a diffusion coefficient D A ) between the locations in an infinite d-dimensional space. We shall denote in the following sections a single A agent as a and a single B agent as b. In between a jumps the B density at the a location grows exponentially (Shnerb et al., 2000a) as B(t) ∼ B(0)e(λ−µ−2d D B )t . The estimation is made by neglecting the flow of B’s from a neighbouring site to the a site. This is justified when the B concentration in the neighbouring sites is much lower than on the a site (it is thus an underestimate of the B growth rate). In the same limit, the ratio between the B density at the a location and the B density on a neighbouring site can be estimated as λ/D B . Consequently, each a jump corresponds to a sudden downwards jump by a factor of λ/D B in the B density. As there are on average 2d D A such jumps per unit time, the net effect of proliferation, diffusion and death, gives the B concentration at the a site as a function of time: B(t) = B(0)e(λ−µ−2d D B −2d D A ln(λ/D B ))t .
(9)
The approximation is in good agreement with the simulation shown in Fig. 3. The slope of the B island maximal density, on a log scale, is indeed ln(λ/D B ), the time dependence of the height of the B island in between a jumps is indeed given approximately by an exponent, (λ − µ − 2d D B )t. Consequently the dashed line (in the inset), which represents equation (3), follows closely the actual growth seen in the simulation (full line). The difference between the theory and simulation is mainly due to the case that two or more a jumps follow each other rapidly, in this case the island’s shape does not stabilize before another a jump is made. These rather rare events, modify somewhat the actual result§ . § The above analysis turns void if D = 0, where the spatial dimensions of the island do not grow B
at all. We do not consider this singular case in this paper. In all other cases, our estimate is an underestimate.
382
Y. Louzoun et al.
Figure 3. Single A approximation. The profile of a B island as a function of time as it follows the random motion of a single A agent. The cross-section of the island is taken through the current location of the A agent. The inset shows the time evolution of the height of the B concentration at the point at which A is currently located (solid line). The B colony grows, although the average growth rate over the entire space is negative (A is extremly low since there is only one A agent in the whole simulation space, thus λA − µ ≈ −µ). The dashed line shows the exponential growth with coefficient λ − µ − 2d D B − 2d D A ln(λD B ). The slope between ‘A jumps’ is λ − µ − 2d D B .
One may ask: what is the situation when B colonies arising around a single a are unstable [i.e., when the exponent in equation (9) is negative]. One possibility is that in such a situation the continuum approximation is valid and the B concentration decays to zero. Another possibility is that, although colonies induced by a single a are unstable, cooperative effects give us back the survival feature. 3.2. Low diffusion approximation. Cooperative effects are clearly observed in the case of a very low A diffusion rate. If the first terms of the exponent are negative (λ − µ − 2d D B < 0), while the last term (−2d D A ln(λ/D B )) is small, a very small nonzero A population would reproduce the exponential growth. For such a low diffusion rate, we can assume that the A configuration is fixed over the time B interactions takes place. For a large enough space, the A will have a Poisson distribution. For a large enough volume, for any average A and for any value m there will be regions in space where Ai >= m (where Ai is the density of A agents at the point i). The m . Thus, if we take probability to have m as at some point is p( Ai = m) = e− A A m! the size of our space to be V e A m! A−m , a large number of points will have at
Proliferation and Competition in Discrete Biological Systems
383
Prediction of MS in the regime of low A diffusion rate
10 12 10 11
B population
10 10 10 9 10 8 10 7 10 6 10 5 10 4
0
10
20
30 Time
40
50
60
Figure 4. The dynamics of the B population in the low A diffusion regime. The parameters used for this simulation are: A = 6, λ = 0.1, µ = 1.0, D B = 0.1, D A = 0.001, λA − µ = −0.4, λ − µ = −0.9. Note that since the A distribution is nearly static the exponent of the B dynamics is constant.
least m as. At every such point we can reproduce the ‘single a approximation’ but with an effective growth rate of λm − µ − 2d D B . If we take m to be big enough then the effective growth rate will be positive and the ‘single a approximation’ will predict an exponential growth. The effective A diffusion rate is m D A , since it is enough for one A to leave the point with the high A density for the effect to fail, but we assume that D A is small enough, and that this effect is negligible. Note that the limiting factor on the value of m is the volume of our space and not only the value of A. Thus even when λ A − µ 0 and λ µ one obtains an exponential growth rate (Fig. 4). 3.3. General theory. We have shown that PDEs fail in two simple to understand cases. The first case represents a low A concentration. In this case the fluctuations are important in comparison with the mean. The second case is a low A diffusion rate. In this case although the A fluctuations are small, they are stable for a long time. We are now in a position to test if we can apply these results to a more general case, where A is not small and D A is the same order of magnitude as the other parameters. We will now show that in a very nonintuitive way the PDE analysis fails even in these cases. This emerging phenomena is completely unpredictable from the microscopic point of view (Fig. 5). In order to understand this phenomena let us decouple the A and B dynamics. If at a given time there is a large A aggregate with m as at a given point, the log of
384
Y. Louzoun et al.
Figure 5. B dynamics within the regime of high A density, mid values of A and B diffusion, λ < µ and λA < µ. The upper plots represents snapshots of the log of the B population at the times pointed by the arrows. One must note the large island that keeps growing through all the snapshots, swallows neighbouring islands and eventually contains a large part of the entire B population, while other islands disappear. The probability that such an island should disappear is inversely proportional to its size. Once a large island is established it has a very low probability of disappearing. However in a limited space after a long enough time even this island will disappear. On the other hand for any given time there is a size of the system for which the total B population keeps growing.
the B population at this point grows following ln(B(x)) = (mλ − µ − 2d D B )t.
(10)
If the local A concentration is not constant, then the B population at this point evolves as: t +t ( A(x, t )λ − µ − 2d D B )dt . (11) ln(B(x)) = t
The log of the ratio between the B population at a small distance r from the point x and the B population at x is: ln(B(x, t)) − ln(B(x + r, t)) =
λr . DB
(12)
Proliferation and Competition in Discrete Biological Systems
385
Let us follow the point of maximal B population in this island. We will denote the original point of maximal B population as x. As long as the integrand on the right-hand side of equation (11) is positive the B population will increase in this site. Once this integrand becomes negative the B population will decrease. However if within the island there are other points where this integrand is positive, the B population can increase in these other points. If the distance of these points from the original B maxima is r, then ln(B) at this point will be at most DλB r lower than in the original position. In order for ln(B) to grow it must increase (on average) in its new position more than DλB r. Condition (12) does not involve the B population. It can be stated entirely in the language of the A dynamics as: t +t λr ( A(x, t )λ − µ − 2d D B )dt > (13) DB t t +t µ + 2d D B r dt > (14) A(x, t ) − ⇒Z= λ DB t t is the interval where the integrand in the left-hand side of equation (14) is positive. Note that the integral itself is defined to be always positive. For the DB , and t is a time interval sake of simplicity we can rescale A to a ∗ = A − µ+2d λ ∗ where a (t) ≥ 0. Equation (14) then becomes: tfinal r a ∗ (t)dt > . (15) DB t0 Let us now analyse the average evolution of an island with a radius of R0 . The average value of r is smaller than R0 . We can thus replace r by R0 in order to obtain a sufficient condition for the island growth. We then use the following algorithm. At every time step we look for the point, within the radius R0 , where Z is maximal. This point can change since the value of a ∗ is changing as the A agents diffuse. Let us name this point the focal point. As long as the value of a ∗ is positive in this point, the B population of the island will grow. Eventually the value of a ∗ at the focal point will become negative. We then reset the value of Z to zero at all the points on the grid, and look again for the point of maximal Z . Note that the point that will become the new focal point can have a negative a ∗ value at the time of the jump, so there can be a small ‘wait period’ between the growth at the previous focal point and the new one. In order for the total B population to grow the average value of Z at the focal points should be bigger then DR0B + wait µλ , where wait is the average wait time. We can thus replace the single A condition equation (9) by a ‘single island condition’: tfinal R0 µ ∗ (16) a (t)dt ≥ + wait . DB λ t0
386
Y. Louzoun et al. 80 70 60
50 40 30 20 10 0
0
0.05
0.1
0.15
0.2
0.25 DA
0.3
0.35
0.4
0.45
0.5
Figure 6. Estimates of Z − wait for a two-dimensional grid with A ∗ = A − 3 and an island diameter of 35. The balance is between the jump needed at the end of the growth period and the growth obtained during the growth period. For any value of D B , one can compute the value of D A for which Z − wait ≥ DR0 . Note that the figure represents B the estimate for a given ‘B island’ diameter, however the estimate obtained from a single diameter is close to the result of the simulation describing the full A and B dynamics.
The new condition is size dependent. In a large enough system, islands of all sizes will temporarily occur. If an island is stable at a given radius it will ‘survive’. In order for the B population to grow, it is thus enough for a single radius to be stable. We check the stability of B islands by testing numerically condition (16) for multiple radiuses¶ . Let us take, for example, a case where λ is half of µ, 2d D B is equal to λ and A is one. In other words the B diffusion rate is as high as the reaction rate, the death rate is much bigger than the birth rate and the A population is not small. According to our scaling relation a ∗ = A − 3. We simulated an island diameter of 35 in two and three dimensions. In two dimensions condition (16) holds for an A diffusion rate of 0.15 and a B diffusion rate of 1.0 (Fig. 6). For three dimensions (not shown) the same condition would hold for A diffusion rates of 0.05. This estimate is obviously an underestimate, since we checked a single diameter, and our algorithm is suboptimal. However simulations with the B agents (Fig. 7) indicates that the computed threshold is reasonable for a system of limited size (100 × 100). ¶ The minimal value of the integral can be computed analytically by ignoring the probability that A
agents shall return to the focal point, and assuming only diffusion away from this point, as we have done in the ‘single A approximation’. However the quality of our estimates depends very much on the return of As to their original position.
Proliferation and Competition in Discrete Biological Systems
387
12
10
10
Total B population
10
8
10
DA=0.1 6
10
4
10
2
DA =0.2
10
0
10
0
50
100
150
Time
Figure 7. Average B cell population on the two sides of the threshold predicted by equation (16). One can see that although the threshold was obtained for a single ‘island’ diameter, it is still a reasonable estimate for systems limited to a size of the order of magnitude of this diameter. When one sets the A diffusion to be lower than the threshold the B population grows indefinitely, while above the threshold the B population decreases to zero.
We have previously shown in the context of a renormalization group analysis that in the regime of high A density and low λ the probability of B survival is much higher in two dimensions than in three dimensions (Shnerb et al., 2000b). This difference is due to the low probability of an a to come back to its original location in three dimensions or more. The current analysis provides an insight on the microscopic details behind our previous results.
4.
M ECHANISMS L IMITING THE P OPULATION G ROWTH
The proliferation of the B agents in the models of Sections 2 and 3 must be limited by some saturation mechanism. We have shown that a high death rate (even a death rate higher than the birth rate) is not enough in order to limit the population proliferation. The required mechanisms can be either an external regulator (Jiang and Chess, 2000), or a limit on the resources available to the proliferating agent (Darwin, 1859). 4.1. Local competition. Many of the mechanisms limiting the B proliferation can be described as competition. Competition for a resource means that whenever two B agents need the resource either to survive or to proliferate only one of them will get it and the other will die (or fail to proliferate) (Fig. 8). In ecology
388
Y. Louzoun et al.
Figure 8. Reaction scheme for a system containing local competition. This system has constant proliferation and death rates as in Fig. 1. However it has an extra mechanism of local competition: Whenever two B agents meet there is a constant probability that one of B. them will be destroyed; B + B
competition can be over food or water, while in the immune system competition can be over access to an antigen. This system is an extension of the logistic equation and can be described using a PDE: ˙ = (λ A(r) − µ)B(r) − γ B(r)2 + D B ∇ 2 B(r). B(r)
(17)
For low B values this system will behave similarly to the system described in the previous sections, where we have already shown the failure of the PDEs. One can analyse this system in two cases. If (λ A − µ) > 0 the total population increases at almost any point in space until limited by the saturation, and then . The discreteness of the A induces stabilizes around a steady state of: B = λ A−µ γ only a low level of fluctuations around this steady state. The more interesting case is when (λ A − µ) < 0. We have shown in the previous sections that, even in this case, the point B = 0 is not a stable steady state (in contradiction with the PDE result). The steady state occurs when the B population in the regions of high B population is saturated by the competition term. The only difference between this case and the one presented in the previous sections is the limited size of the high B population zones (Fig. 9). Paradoxically the conditions where PDE and the MS show similar features are seldom realized in nature. They represent extreme cases where a certain agent is able to fill all the available space. In the animal world only a few species enjoy such a status (for example, humans and their parasites). In immunology this regime corresponds to a cell type filling all the blood and immune systems (systemic diseases). Except for these extreme cases the situation where (λ A − µ) < 0 seems to be more appropriate for realistic systems.
Proliferation and Competition in Discrete Biological Systems
389
Figure 9. Comparison between PDE and MS for the system defined by Fig. 8. This system has an average negative growth rate yielding an exponential decay in the PDE description as in the previous section. In the MS we observe a limited growth of the B islands, which is due to the competition between Bs. Thus even regions where the local B concentration would grow are limited to a finite concentration.
4.2. Global competition. The proliferation of the B agent, and the competition between different B agents can occur over very different spatial scales. The proliferation is a local mechanism; the factor leading to competition can be global. In ecological systems the global character of the limiting factor may be related with the fact that the predators move much faster than their prey and that their population is proportional to the total prey population in their hunting region. Therefore, their population constitutes an indirect interaction between distant prey locations. We will show that in such cases the inhomogeneous effects due to discreteness are even more important than in the previous sections. As an aside let us note that this kind of reactions although formally local (a lion can eat a prey only if they meet) is very effective in inducing top down effects; that is the regulation of the local prey population is enforced by their total population. MS can describe in a natural way only local interactions. It is nontrivial to describe with MS factors acting on widely different scales. The solution in this case would be to create a new type of Hybrid Model, which will allow in particular to express the interaction between agents and their averages over various ranges.
390
Y. Louzoun et al.
Figure 10. Global competition reaction scheme. We added to the reaction described in Fig. 1 a global competition interaction. The death probability of each B agent is proportional to the total B agent population.
A simple example for such a system is obtained by modifying the regulation mechanism of the system presented in the previous section (Fig. 10). There the two mechanisms proliferation and competition were local interactions. We can now replace the local competition by a global competition. Even though from the microscopic point of view the nonlocality might look unnatural this is the only way that nonlinearities can arise at the macroscopic level. Indeed an ODE term of the form B2 can appear only as the average over terms of the type B(r) ∗ B. Therefore the modification of equation (6) would be ˙ = λ A(r)B(r) − µB(r) − γ B(r) B + D∇ 2 B(r). B(r)
(18)
That is we replaced the local competition in the system described in Section 4.1 by a global competition term. We find that passing from local to global reactions has crucial implications: by comparing the behaviour of the system equation (18) with the system equation (17) one finds that they are dramatically different. Even in the region where the system equation (17) is described correctly by a PDE, a PDE is still incapable of describing properly the system equation (18). Indeed in system equation (18) (for most parameter sets) the entire space is empty of Bs except for a single island centreed around the point where the A concentration is maximal (Fig. 11). We can understand this system by ignoring (temporarily) the diffusion term and simplifying the equation to : B˙ = A(r)B(r) − B(r) B = B(r)(A(r) − B).
(19)
We rescaled the equation to λ = 1 and included the death term µ into the A density. This can be done with no loss of generality.
Proliferation and Competition in Discrete Biological Systems
391
Figure 11. Comparison between PDE and MS for a system with an average negative growth rate and a global competition term (as defined in Fig. 10). In this model a single large island with a very high B concentration appears. The competition with this island destroys all the other small B islands. The (a) and (b) drawings show the log of the B agents concentration at different times. The (c) drawing shows the comparison between the PDE and the MS.
One can see that in every point that A(r) is lower than Amax the B(r) density will vanish, while around the point rmax defined by A(r) = Amax a region of high B density will appear. Suppose that one starts with all the B(r)s very small; as long as the average concentration of the Bs ( B) is lower than Amax , the local concentration B(r) at the point rmax will rise [according to equation (19)]. This in turn will raise the value of B. Eventually B will reach the value of Amax . What will be the equilibrium value of the B(r)s in every other point? In every point where A(r) < Amax , one will have B(r)( A(r) − B) = B(r)( A(r) − Amax ) < 0, which according to (19) means that the concentration at those points will decay to 0. It is clear therefore that the entire B population responsible for B = Amax is due to a very dense island situated around rmax . The equality B = Amax is not affected by the B diffusion, which only has the effect of spreading the distribution of the B(r)’s in a limited region around rmax . One sees that the dynamics is dominated in this case by the most extreme stochastic local fluctuation in the A(r) distribution, and not by the global or average properties of the A(r)s.
392
Y. Louzoun et al.
Figure 12. The model defined in Fig. 9 in a regime of high A diffusion rate. When the A diffusion rate is very high, a few large B islands are created. As oppossed to Fig. 11, the largest B island does not always dominate, and smaller island temporarily appear and disappear. The various islands emerge an adaptive behaviour: they look for the maxima of the A concentration, split, merge and die. The top windows shows different snapshots of the B distribution. The bottom window is the long-term evolution of this system showing intermittent fluctuations.
In the presence of A diffusion two regimes appear. In most of phase space all the B population is concentrated in a single island, but in a small part of phase space, multiple islands can exist for a long time (Fig. 12). The phase transition occurs when the destruction time of small B islands is longer than the time over which the λA −µ−2d D B , largest island is stable. The total B density is determined by Btot = big isl γ where Abig isl is the A density at the centre of the largest island. The disappearing 1 1 = λ( Abig isl −A . time of smaller island is proportional to λ Asmall isl −µ−2d D B −γ Btot small isl ) The difference between the A density in the largest islands and the next smaller islands is of the order of unities. The largest island will thus lose its advantage over other islands at a rate proportional to D A . Multiple islands can thus exist if DA ∼ DλA > 1, on the other hand if D A is much bigger than λ then λ( Abig isl −Asmall isl ) no island will survive and the results will converge to the PDE estimate. We are thus left with a small part of shape space close to the phase transition in which multiple islands can survive. Within this regime, the diffusion of the A does not
Proliferation and Competition in Discrete Biological Systems
393
lead to smoothing of the B fluctuations. It only leads to a movement of the high B population zone. This zone does not disappear, since its lifespan is much larger than the lifespan of a single B (it is the time it takes to destroy a macroscopically large number of Bs). Within this lifespan a new large fluctuation of A will appear within this zone, and this zone will start to grow again. The b islands population is then greatly affected by the local changes in the a distribution and the entire system undergoes intermittent macroscopic fluctuations [in certain systems such fluctuations can have catastrophic consequences (Louzoun et al., 2002)]. 4.3. Emergence of complex behaviour. In the first sections it has been shown that the discrete character of the elementary agents A leads to unexpectedly rich and complex behaviour of the B population. The complex spatio-temporal emerging structure of the B population cannot be predicted by a PDE approach. In the present chapter we introduced a few types of competition terms. Local competition limits the size and maximum concentration of each macroscopic B island. Rather than interfering with the emergence of macroscopic features, the competition terms turned out to enhance the local and complex character of the B islands. Global competition can endow these macroscopic islands with an emerging social (more precisely anti-social) behaviour. The evolving macroscopic objects now have apparent ‘goals’: seeking food (As), multiplying, dying or destroying the competing islands. The lifespan of the emerging islands is many orders of magnitude longer than the lifespan of their components (which of course remain nonadaptive and totally mechanical during the entire history of the system).
5.
D ISCUSSION
Dynamic biological systems were traditionally described using PDEs. PDEs were invented and proved useful in simple physical systems that fulfilled certain restrictions: • A very large density of interacting particles. • A diffusion rate high enough to smooth local fluctuations. • An immediate result of each interaction, i.e., no delay between the interaction and its result. These assumptions are often false in cellular systems. • The density of interacting particles is often small; although the total number of interacting agents may be very large, the local number of interacting cells can be very low. In order to use a PDE approximation the Poissonian fluctuations in the number of interacting agents should be much smaller than the average number of interacting agents. Especially in the case of autocatalytic systems the density should be high enough that the probability of a large
394
Y. Louzoun et al.
fluctuation should be low. Note that in contradiction with the naive intuition, the probability of larger fluctuations increases with the size of the modelled space. • The diffusion rate can in many cases be low. The diffusion rate should be defined as the ratio between the rate of the random motion of the agents and their interaction radius. Moreover, compartmentalization such as physical objects in ecology or endoplasmic reticulum in intracellular interactions can further limit the diffusion. To summarize, PDEs can be used only if the agents move freely faster than the interaction radius divided by the reaction rate. • Most biological processes require the production of new molecules or the division of cells. Each of these processes can take between hours and days. Thus, there is a delay between an action and its result. This aspect is the most problematic since even minor delays can change the dynamics of the system. A simple example is the loss of stability of the fix point in predator– prey in the presence of a small delay. At the current point we cannot provide a ‘thumb rule’ for when delays can be ignored. We studied in this paper one of the main sources of complexity: the combination of microscopic inhomogeneity and autocatalysis (Shnerb et al., 2000a). In biological systems, there is an inherent source of microscopic inhomogeneity—the discreteness of the entities composing biological systems. Moreover, most biological systems involve proliferation, which is the basic autocatalytic mechanism. We showed that in some very simple cases containing these two elements PDEs fail to describe appropriately the dynamics. This is as well: systems, which were too simple to present complex emergent features in the framework of PDEs turned out to present when studied appropriately a complex behaviour. We used in addition to the straightforward MS intermediate techniques that allow to take explicitly into account the intermediate relevant dynamical scales between the MS and PDEs. In the systems simulated in our paper, autocatalysis results in the flow of information from small-scale to large-scale features; from the microscopic details of the system to its macroscopic behaviour. The rapid evolution of small-scale features to a macroscopic scale affects the distribution of the macroscopic variables. In particular the macroscopic variables which are dominated by the sum of a small number of large values (instead of the sum of a large number of small values) have a very specific (fractal-power law) space–time behaviour. Such effects are enhanced when there is a direct effect of the macroscopic values on the microscopic mechanisms. Many biological systems evolve at time scales much longer than their fundamental microscopic time scales. In PDEs the emergence of long time scales cannot occur naturally, while in MS they emerge spontaneously. In MS, the agents selforganize (Atlan, 1987) in collective structures with emergent adaptive properties and with dynamical space–time scales of their own. The effective dynamics of these emerging objects is much slower than the one of their components [and this fact constitutes their very definition (Solomon, 1995)]. Analysing the system at the
Proliferation and Competition in Discrete Biological Systems
395
level of the emerging objects automatically requires the transition to longer time scales and larger spatial scale. Each level corresponds to a new space–time scale in a hierarchy of complementary representations of the same complex macroscopic system. To summarize, we found in this paper that one of the core mechanisms for the emergence of collective macroscopic objects is the interplay between autocatalysis and fluctuations inherent to the discrete structure of the microscopic components of the system. PDEs or even stochastic PDEs cannot describe this mechanism, since both the emergence and the disappearance of these collective objects are a result of their discrete nature.
R EFERENCES Atlan, H. (1987). Self-creation of meaning. Phys. Scr. 36, 563–580. Boyce, W. E. and R. C. DiPrima (1986). Elementary Differential Equations and Boundary Value Problems, 4th edn, New York: Wiley. Cohen, I. R. (2000). Tending Adam’s Garden: Evolving the Cognitive Immune Self, San Diego, CA: Academic Press. Croft, M. (1991). Activation of naive, memory and effector T cells. Curr. Opin. Immunol. 6, 431–437. Darwin, C. (1859). The Origin of Species, London, U.K.: John Murray. Jiang, H. and L. Chess (2000). The specific regulation of immune responses by CD8+ T cells restricted by the MHC class Ib molecule, Qa-1. Annu. Rev. Immunol. 18, 185–216. Kauffman, A. S. (1995). At Home in the Universe: the Search of the Laws of Complexity, U.K.: Oxford University Press. Levy, M., H. Levy and S. Solomon (2000). Microscopic Simulation of Financial Markets; From Investor Behavior to Market Phenomena, New York: Academic Press. Louzoun, Y., D. Mazurski, J. Goldberg and S. Solomon (2002). The risk of being unfair, the unstable fate of globalisation (submitted). Louzoun, Y., S. Solomon, H. Atlan and I. R. Cohen (2001). The emergence of spatial complexity in the immune system. Physica A 297, 242–252. Malthus, T. R. (1798). An Essay on the Principle of Population printed for J Johnsonin, St Paul’s Churchyard, London (reprinted by Macmillan et Co 1894). Meier-Schellersheim, M. and G. Mack (1999). SIMMUNE, a tool for simulating and analyzing immune system behavior DESY-99-034. Papoulis, A. (1984). Probability, Random Variables, and Stochastic Processes, 2nd edn, New York: McGraw-Hill, pp. 144–145. Segel A. L. (1991). Biological Kinetics, Cambridge Studies in Mathematical Biology 12, Cambridge, U.K. Seiden, P. E. and F. Celada (1992). A model for simulating cognate recognition and response in the immune system. J. Theor. Biol. 158, 329. Shnerb, N., Y. Louzoun, E. Bettelheim and S. Solomon (2000a). The importance of being discrete—life always wins on the surface. Proc. Natl Acad. Sci. 97, 10322. Shnerb, N., E. Bettelheim, Y. Louzoun, O. Agam and S. Solomon (2000b). Adaptation of autocatalytic fluctuations to diffusive noise. Phys. Rev. E 63, 021103.
396
Y. Louzoun et al.
Solomon, S. (1995). The microscopic representation of complex macroscopic phenomena: critical slowing down—a blessing in disguise, in Annual Reviews of Computational Physics II, D. Stauffer (Ed.), World Scientific, pp. 466. Stauffer, D. (1999). Anual Reviews of Computational Physics VII, Singapore: World Scientific Publishing Company. Verhuulst P. F. (1838). Notice sur la Loi que la Population Suit dans son Accroissement; Correspondence Mathematique et Physique, A. Quetelet (Ed.) pp. 113–121; (English translation in Smith, D. and N. Keyfiz (1977). Mathematical Demography, New York: Springer, pp. 333–337). Received 7 May 2002 and accepted 27 November 2002