The remarkable behavior of a map displaying derived from a simple model for a biological conflict

One of the simplest yet profound mathematical models for biological growth emerged sometime in the middle of the 1800s due to the work of Verhulst. It describes population growth thus: let be the population of the organism at time . is the Malthusian parameter or the rate of maximum population growth. is the carrying capacity or the maximum sustainable population in a given environment. Then we have the following ordinary differential equation (ODE) as a descriptor of population growth:

This is the logistic model that profoundly informs us about various aspects of biology. By define a new we get:

This simplified form can be next converted into a discrete equation thus:

This is the famous logistic map. Simple as it looks, starting with von Neumann’s work and then that of Ulam and Feigenbaum thereafter, it became clear that it exhibits notoriously complex behavior that have considerable implications for dynamics of systems. We had earlier described an exploration of some of this via the vehicle of Kaneko’s coupled-map lattices. Since our teenage years we have been fascinated by the complexity in simplicity of this map and wondered if we might uncover other such cases in simple models of biological systems. Here we describe an example of such, which we discovered.

Consider a two species-system where one species is a parasite or predator, whose growth rate depends on the number of the host/prey species that it utilizes. The host/prey species is communal in that it produces certain public goods, which all members of its community share. It incurs a certain cost for producing these. More its population less is this cost because it gets spreadout over more individuals. It grows like any other free-living organism otherwise with growth rate proportional to the number of individuals at a given instant. Its growth is further negatively affected by the predator/parasite which kills all infected individuals. In principle, this system might be described by the below ODEs:

By simplifying this by taking and expressing it as a discrete equation we get:

Since is a cost of public goods production by the communal organism it will always be negative, while being a regular growth constant will be positive. This is the form of the map, which we shall explore further for its dynamics. A few things are readily apparent if then the first term of second equation of the map goes to resulting in a singularity. Second both and can wander away towards ; hence in out experiments we terminate our runs if certain arbitrarily defined limit is reached.

Figure 1

Figure 1 shows a run of this system for 100 generations with and initiated with . The expected pattern of tracking is seen in the oscillatory pattern typical of predator/parasite-prey/host system. We can better understand this by plotting a phase diagram of as can be seen in Figure 2.

Figure 2

In Figure 2 we have 9 different runs of the map with different starting and evolving up to 10000 iterations or terminated if it hits a singularity (in practice we set some absolute value cut off of say 100 above which we stop further iterations). The maps are symmetric about the line, which arises from the form of the first equation in the map. The map is always bounded by an ellipse of the form , where . Hence, we call it the 2-parameter elliptical map. Figure 2 reveals that the map is extremely sensitive to changes in the initial conditions: changes in the 3rd to 5th place after the decimal point can result in dramatic differences in the evolution of the .

Figure 3

Figure 3 shows 9 runs of the evolution of for different . We gain see that the map is very sensitive to changes in both the parameters and the evolution of the initial point dramatically varies with these changes.

Figure 4

In Figure 4 we run the following experiment: We chose a square region defined by the diagonal (0.2, 0.2):(-0.2, -0.2). Within this square we chose 3 different sets of 100 random points and allow each set of these points to evolve for a given . Thus, each row in Figure 4 is the evolution of the three different sets of random points for a given . The evolution of each of the 100 random points for a given run is plotted in a different color. We observe that each run shows a different color pattern because the points are random and evolve very differently, as we know from the single point runs shown in Figures 2 and 3. However, remarkably, the overall structure of the map for a given converges to a constant form. This suggests that even though the map is very sensitive to the initial conditions, the average evolution of a given region, i.e. a random set of points drawn from the region, is conservative for a given parameter set.

Figure 5

In the next experiment (Figure 5) we try to understand the effect of the on the average evolution of a given “region”. However, we use a slightly variant definition of the region in this case: we chose 64 points on a circle of radius each separated from the next by an angle of radians from 0 to . We then let these points evolve for 5000 iterations or until they explode to a singularity. The evolution of each point is plotted in a different color as above. For each row is constant while changes. For each column is constant while changes. The effect of increasing is like that of a magnifier. Thus, the increasing acts like zooming in on the structure of the interior of the map. The obvious effect of changing is immediately apparent: with increasing the eccentricity of the bounding ellipse of the map increases. As the bounding ellipse is of the form , essentially . also defines the internal structure of the map, such as its basic geometry as well as the zones of complete occupancy and the zones of exclusion. One notices these changing with changes in .

Figure 6A. 0.61, , .62, .99, 1, 1.01, 1.24, , 1.2475

This effect of the role of on the internal structure is explored at greater depth in the several parts of Figure 6, which were essentially rendered by running the map as in Figure 5 with a constant and various . Remarkably, at , where is the Golden ratio, the map assumes a 10-rayed form with a tendency to diverge rapidly only along those rays (panel 2 of Figure 6A). This value is the ratio of the side of a regular pentagon to its diagonal. This is the central value of the “pentad” state, i.e. close to it on either side will show the characteristic pentad geometry as can be seen in panels 1 and 3 of Figure 6A. Similarly, we observe that is the central value of the hexad state (panel 5 of Figure 6A). At values close to it on either side the map assumes a hexad form as can be seen in panels 4 and 5 of Figure 6A. Another central value is . This number is derived from heptagonal equivalents of the Golden ratio: it is the ratio of the length of the larger diagonal to the smaller diagonal of a regular heptagon. Here again we see a radiating structure with 14 rays (panel 8 of Figure 6A). As approaches this value we see the interior of the map assuming a heptad structure (panel 7 and 9, Figure 6A). These cases illustrate that for values close to the central value but less than it the map displays a -ad polygonal interior structure. For the values close to but greater than the central value the -ad interior structure assumes a rosette form with zones of exclusion defined by the rosette.

Figure 6B. where 1.4, , 1.53, , 1.61, , , 1.73,

Continuing this way, we can find further central values associated with higher -ads. At , the map displays 8 rays (panel 2, Figure 6B). This corresponds to the ratio of the largest to the shortest diagonal of a regular octagon. As we approach this value we get a map with an octad structure (panel 1, Figure 6B). This is followed by the central point defined by , which displays a map with 18 rays (panel 4, Figure 6B). This corresponds to the ratio of the largest to the shortest diagonal of a regular nonagon. Thus, as we approach this value the map takes on a nonad structure (panel 3, Figure 6B). The next central value is (panel 6, Figure 6B), which again results in a 10-rayed map and values approaching it have a decad structure (panel 5, Figure 6B). The hendecad (11-fold) central value is , which results in a 22-rayed map (panel 7, Figure 6B). Similarly, is the central point for the dodecad structure (panels 8 and 9, Figure 6B). Studying these central values of , which mark the -ad geometry of the map, we can derive a general relationship between these to the ratio of the diagonals of the corresponding regular -sided polygon thus: Given a vertex of the regular polygon with sides, let be the
first side of the polygon emanating from that vertex. Then will be the successive diagonals and will the second side emanating from that vertex. Then,

Using some trigonometry, we can show that for a -sided polygon with unit sides the length of the th diagonal (including sides) connected to a given vertex is:

This defines the values at which the 2-parameter elliptical map will take a divergent form with -rays (if is even) or -rays (if is odd). At values approaching these it will show a pronounced -ad internal structure. This explains why the pentad is centered at where and for a unit pentagon. The tetrad has no ; hence its central value would correspond to : that’s why we see the tetrad structure at closer to zero (Figure 5). For we already reach . As increases we have:

Thus, the higher -ad structures are closely clustered near 2. The derivation of the above formula for enabled us to investigate other central points, where instead of integer we have a rational number ; being mutually prime. Thus we consider of the form:

Figure 7A. 0.5165, , 0.5195, 0.8672, , 0.8683, 0.39, , 0.3911

The transition points associated with such rational values are explored in the various parts of Figure 7. In Figure 7A, we show the maps associated with (panels 2, 5, 8). Each of these shows the -rayed structure typical of the the central points. Further, the number of rays corresponds to the value of . Again as in the case of the corresponding to integer as we approach the corresponding to we get -ad structures with a similar pattern as with the former. Here we have chosen 24, 28 and 32 (the syllable counts of the gāyatrī, uṣnik and anuṣṭubh meters).

Figure 7B. , -0.44, , 0.342, , 0.24, -0.28, , -0.2855

In figure 7B we explore further cases of corresponding to different . Notably, using these we get lower values of at which we can get central points of hept-, hendec- and 13- ad structures. Interestingly, at some of these values is negative (Figure 7B, panels 1,2 and 7,8,9). When is negative, we see that the map becomes a mirror image of the positive , now being symmetric along the line and maintaining an elliptical form between . A negative is not quite meaningful in the original biological model in which the map was defined because it is a growth constant and positive. However, we could still contrive a definition wherein the negative specifies a decaying population, which has a positive effect by reducing the parasite load, which is directly dependent on it, and thereby allowing relative growth.

Figure 8. , , ,

In Figure 8 we explore the more complex effects of the fraction on the structure of the map. The outer structure of each of these maps respectively shows 64, 128, 96 and 100 rays as would be expected from value of . However, in the case of , a closer look at the interior of the map reveals first a 17-ad and then closest to the center a tetrad structure (Figure 8, panel 1). Similarly, with we see an inner 33-ad and innermost tetrad structure (Figure 8, panel 2). For we likewise see a inward progression first to a 25-ad and then a tetrad structure (Figure 8, panel 3). For we see an inward progression to a 13-ad then tetrad structure. In each of these cases the innermost state is a tetrad. The corresponding value of is: 0.1960, 0.0981 0.1308, 0.2506. In each case is much smaller than and closer to 0, which is the center of the tetrad state. Hence, it in each of these cases at the innermost level the effect of tetrad central value i.e. 0 dominates. Beyond that we see that is close to , thus resulting in the emergence of 17-ad structure in next level. Similarly for , we see that it is close to which results in the 33-ad structure at the second level. It is notable that both have a comparable second layer -ad of the form . Now in the case of we have close to it resulting a 25-ad second layer. Finally, with we have near it resulting in a 13-ad second layer. Thus, the value of in the fraction in the formula determines the general -ad structure of the map and actual value of determined by the fraction determines its fine structure. In general terms, the resultant map at a certain may be viewed as a superposition of the corresponding -ad structure in the outer layers with all the -s of the nearby fractions competing to contribute to the inner layers. Thus, away from the strong central values, such as those from integer we get more complex maps.

The lower values at which we get higher -ad structures result in lower eccentricities of the bounding ellipse of the map. This allows us to visualize them better than at values derived from integer because the latter result in bounding ellipses with high eccentricity. The absolute value of determines the eccentricity of the bounding ellipse of the map: As (the eccentricity of an ellipse ; at eccentricity 0 we get a circle and at eccentricity of 1 a parabola). We can derive the below formula giving the relationship between the the parameter and . For the bounding ellipse with corresponding to the discussed above we can write the eccentricity thus:

By expressing in terms of using the above formula we derived for it and doing some trigonometric manipulations we get:

Thus, one can see how only for , and the map is elliptical.

In a final experiment we performed a Julia-like operation: was kept constant and was changed. The evolution was studied for all points separated from each other by 0.001 in the square defined by the diagonal (-.2,-.2):(.2,.2). The points were allowed to evolve to a maximum of 5000 iterations if they did not crash to a singularity or reach an absolute value greater than (100,100). Each point in this square is colored by the number of iterations it survived on a color scale from black to gray via the colors dark blue, dark red, burlywood, coral, chartreuse, and cadet blue. Thus, the points that crash out after the first iteration are black and those that last all 5000 iterations are gray. Figure 9 shows the results of this run for selected

Figure 9. , .25, , .26, .6, , 0.62, 0.99, , .62, .99, , 1.05.

From this figure it becomes clear that beyond the basic symmetry axis along the line the plots are rather “entropic”, i.e. there is not much sign of order in arrangement of the points in terms of the number of iterations they survive. However, one tendency becomes clear, consistent with the actual structure of the maps plotted above: If takes the form with integer , like (Figure 9, panel 5, 8) then the majority of points do not survive the entire 5000 iterations — we see little gray. Thus, when corresponds to or fractions of the form with relatively small then most points are not stable in their evolution and crash to a singularity or diverge rapidly. On the other hand, if is defined by a with relatively large (Figure 9, panel 2) or is not defined by a rational number in the formula, then we get many more points evolving stably over a large number of iterations (the remaining panels of Figure 9) — here we see much more gray. The fraction of points surviving at least 5000 iterations is approximately half the total number of points (fraction ranges from 0.509 for Figure 8, panel 7 to 0.599 for panel 2). The remaining points bail out at fewer iterations in decaying fractions when grouped by bins of 100 from 100 to 4900 iterations (Figure 10). However,
the mostly disordered distribution of these points, relates to the mostly convergent structure of the map for a given region that was illustrated above.

Figure 10. Histogram of bailout or survival till 5000 iterations with same parameters as in Figure 9.

In conclusion, while this simple map shows astounding visual beauty, interesting mathematical structure, and deep links to polygons and coprimality of numbers, one may ask if it provides any biological insights at all. At first sight it does seem like a rather unrealistic model for a biological system given the negative values of population that it generates. But perhaps even this can be remedied by the suggestion that these negative values are not absolute population sizes but only represent the part of a larger population which actually participates in the dynamics. But the key lessons the map provides for biological systems are of a more qualitative type: 1) On one hand it shows how systems even with just two species locked in biological conflict can show chaotic dynamics, where the long term out come is highly dependent on the initial conditions. Even small changes can push the system, for the same parameter set, from long-term stability to instability (defined as singularity or blowout divergence). 2) On the other hand, despite the above, it shows that, as an aggregate, a set of different initial conditions strongly tend to converge to dynamics of a constant form for a set of parameters. Thus, even if individual starting populations have very different fates, a population of populations tends to have predictable dynamics when taken as an aggregate. Such analogies might help understand the stable conflict despite ongoing arms races that have lasted for few billion years (e.g. between bacteria and their phages). 3) The high instability of certain points such the corresponding to integer and certain rational fractions suggests that the dynamics of the system might have some deep numerical constraints and natural selection might essentially have to play the role of selecting between parameters that conform to underlying numerical constraints that favor long terms stability.