Nonlinear self-organized population dynamics induced by external selective nonlocal processes.

Self-organization evolution of a population is studied considering generalized reaction-diffusion equations. We proposed a model based on non-local operators that has several of the equations traditionally used in research on population dynamics as particular cases. Then, employing a relatively simple functional form of the non-local kernel, we determined the conditions under which the analyzed population develops spatial patterns, as well as their main characteristics. Finally, we established a relationship between the developed model and real systems by making simulations of bacterial populations subjected to non-homogeneous lighting conditions. Our proposal reproduces some of the experimental results that other approaches considered previously had not been able to obtain.

Introduction

Broadly speaking, a population represents a group of individuals sharing the same space, being possible to find, for instance, groups of viruses, bacteria, or people. The evolution of these groups, especially the human population, has fascinated scientists for centuries. Questions related to the sustainability of our species led to developing models such as the exponential growth of Malthus [1], or the logistic growth of Verhulst [2], which sets a limit to the population size, the so-called carrying capacity. Although they are not entirely realistic, these models are significant, since in some stages of human development, the population size has presented a marked exponential character, being practically stationary in others [3], [4]. So it could be hard to make an accurate prediction about populations’ behavior, requiring the constant improvement of the models used.

The influence of germ populations on people’s lives is particularly evident. With some frequency, new viruses arise for which we have no defense, spreading around the world at an incredible speed, the most recent case being the coronavirus SARS-CoV-2, which gives rise to the disease known as COVID-19 [5], which has radically changed our way of life. Similarly, bacterial populations often exhibit mutations and collective behaviors that can produce resistance to antibiotics [6].

Examples such as those described explain why there is a great interest in the study of population dynamics. This theme has traditionally been addressed from the establishment of conservation equations that take into account the factors that affect the size of the population studied. Thus we arrive at the reaction-diffusion equations (RDE), which describe the temporal evolution of the population density in the form [3] These equations frequently appear in the study of systems of the most diverse nature [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], leading to interesting phenomena, such as critical behavior, multiple stationary states, spatial patterns, wave fronts and oscillations. Here Γ is the source term, whose specific form will depend on the problem under study. It can create (destroy) individuals in the region Ω. Meanwhile, the second term of the right-hand side of (1) regulates the flux of individuals, being D the diffusion coefficient.

In the present work, we will use the RDE to study the phenomenon of pattern formation. Our proposal, in addition to incorporating known cases, also takes into account aspects such as the non-local character of the interaction between individuals, thus allowing us to reproduce experimental features not previously obtained employing other approaches. The structure is as follows. Section 2 analyzes the model in-depth, leading to results that will be tested later in Section 3, which shows the details of the computational implementation. Section 4 deals with the study of bacterial populations, revealing the physical meaning of some parameters. Finally, Section 5 shows a summary of the main results obtained, as well as some ideas for future research.

A general equation for population dynamics

We propose a set of equations of the formwhere is a functional operator acting on u through the kernel g, being . Here a and b are positive constants representing the growth and competition rates, while α and β are growth and competition length parameters respectively, which characterize the range of interactions in the system. When Eq. (2) has a steady-state () with densityPrecisely this state in which the density is the same for all points of the space will be the state without pattern formation (SWPF). In other words, we will interpret uniformity in density as the absence of patterns. Later on, we will find other states belonging to this category.

Eq. (2) is very general, however, it is not useful unless we define the functional operator . A good and practical operator should meet two main conditions. First, it should be linear, which is more simple, allowing, under certain conditions, the obtaining of a Malthus growth term, as well as a competitive term proportional to u 2(x, t). Second, it should be able to go beyond already known results, exhibiting new features. Thus we propose a functional operator as (from now on, we will consider the one-dimensional case)A relation similar to Eq. (4) has already been used in Refs. [49] and [50] for studying pattern formation.

Substituting (4) into (2) leads towhere and are non-local correlations for growth and competition that weigh the interaction between individuals. Here the nonlocal kernels play a role similar to that of the memory kernel in the Generalized Langevin Equation [51], [52], [53], [54], [55], [56], [57], [58]. In the case of Langevin dynamics, the memory kernel reflects the fact that events occurred in the past can affect the present. On the other hand, due to processes like diffusion, the density on a certain point should be correlated with its neighborhood [13], [16], [49]. We will assume that the nonlocal kernels are even functions normalized to the unit in the domain Ω.

The integrals in (5) can be expanded in series as [49] where . Introducing the first term of (6) into (5) yields the Verhulst equation. The second term () is the diffusive part, while the higher-order ones are dispersive.

Now we can obtain various formulations for the growth dynamics. For example, if we develop the first integral of (5) until the second term, and the second integral until the first one, it is obtained the famous Fisher equation [59]:where the diffusion coefficient is [49], [50] On the other hand, by keeping the second integral in Eq. (5) and expanding the first integral to the second-order, we get the generalized Fisher-Kolmogorov equation [13], [14], [16], [17], [60]. These reaction-diffusion equations are particular cases of (5), which is a particular case of (2).

A simple choice for the non-local kernel isthat allows writing (5) aswhich was solved numerically in Ref. [49], considering periodic boundary conditions (PBC) [61] as well as the initial condition u(x, 0).

Before continuing to work, it is convenient to express (10) in terms of dimensionless quantities [3], what can be done through the following transformations: where L is the size of the system, while the hat indicates that the magnitude is dimensionless. Now (10) takes the formHere all the quantities involved are dimensionless. The hat was discarded for the sake of simplicity in notation. Notice that (18) does not depend explicitly on neither a nor b. At the same time, x, α and β are fractions of therefore, their values will always be in the range [0, 1]. In this way, we can work with (18) and subsequently obtain the results corresponding to systems of specific dimensions through the transformations (12)–(17).

At first glance, it seems that the parameters α and β can assume any value within the range (0,1). However, the fact of using PBC will impose significant restrictions on them. Let’s analyze, for example, the case of a population located in a one-dimensional space of size .

If the value of the parameter that characterizes the range of interaction between individuals is then an individual A 1 located in the middle of the space would interact with everybody. In the case of a member A 2, at the right of A 1, in the position x 2, its interaction interval will be while part of its reach extends beyond the border of the system, specifically the range . But this region, according to the PBC, corresponds with the band . With that, we obtain that the member located in A 2 can also interact with all the population. We get the same result when analyzing points to the left of A 1.

Doing a similar analysis for κ > L/2, we observe that part of the interaction is double-counted. So when using PBC, the values of α and β are restricted to the rangeIn the following paragraphs, we will see that there are restrictions that further limit the values of these parameters, especially those of α.

The condition (19) is the result of applying PBC. With other boundary conditions, the ranges of variation of the parameters characterizing the interactions would be different [62], [63]. The use of PBC in the present work is because they are a common choice in the literature that studies phenomena such as pattern formation [13], [49], [50], [63], [64], [65], allowing quite a good description of experimental results. Calculations in Section 4 use different boundary conditions, which would give us a better idea about their influence on final results.

Stability

Having established an explicit form of the model in (18), it is very important to make the linear stability analysis of the steady-state [16], [49], [66], [67], [68], since it will allow us to determine the conditions for the appearance of patterns. For that we are going to use the functionwhere u 0 is the homogeneous state given in (3), while ξ ≪ u 0 is a perturbation which, for large t, can either grow or die out, depending on the value of η, in such a way thatThis approach has been widely used in nonlinear dynamics [16], [27], [44], [45], [49], [50], [69], [70], [71], [72]. Now we can use the relations (12)–(17), together with to obtain the dimensionless version of (20), namely,where, once again, the hat is discarded to simplify the notation.

Substituting (25) into (18), doing some calculations, and retaining the terms until first-order in ξ, one obtains thatSo we can say that patterns appear whenbeingThe result given in (26) was already obtained in Ref. [49], where it was also reported that, when η > 0, patterns appear if α < β. Meanwhile, the allowed values of the wavenumber k will depend on the boundary conditions, as well as on the form of the solution (see next paragraphs).

Analytical solution for the non-local equation

Now let’s search for the solution of (18) in the formwhere B 0(t) is the contribution to the homogeneous part of the density, while the coefficients B(t), with n ≥ 1, are associated with the inhomogeneous spatial behavior. The periodic boundary conditions allow us to say that which, in combination with (29), leads to . Substituting (29) into (18) we get the following system of differential equations: where being for n odd, and if m ≥ n/2.

Using (29) we obtain that the size of the population at any given time will besince the contribution of the periodic functions is null. So that the size of the population is directly proportional to the value of B 0(t). Nonetheless, B 0(t) depends on all the other coefficients, as shown by the Eq. (30).

In the steady-state (t), Eqs. (30) and (31) are equal to zero, which allows calculating the coefficients as In the absence of patterns, which in terms of dimensionless quantities becomes . This implies that . Considering the linear stability analysis previously performed, we expect the presence of patterns to maximize the population size when the steady-state is reached. So, to observe patterns, it is necessary that 1 < B 0(t) < ∞, which suppresses the negative solution in (35), getting that

In this way, the series is (at least conditionally) convergent. So we can always find a real number μ > 1 so that the general term of the convergent seriesis proportional to the general term of that is,with 0 < λ < ∞. Therefore, the asymptotic form of the coefficients B(t) could be given by

Based on (40), let’s assume that the coefficients are well behaved, in the sense thatand that their absolute values decrease rapidly with the increase of n, so that the sign of is determined by the first nonzero term in the series.

As the aforementioned allows us to establish that the first non-null term of must be negative. Meanwhile, the condition (19) yields f(k 1 β) > 0, leading to which implies that all terms of the sum in (43) are identically null. Since the bracket factor is, in general, different from zero, we conclude that at least one of the coefficients, either B(t) or should be null.

Analyzing the behavior of ψ(t) for different values of n, it is possible to show (see the Appendix A) that, when the index of the first non-null term of the series presented in the Eq. (37) is M, the only non-zero coefficients will be those whose index is a multiple of M, obtaining thatAs the first term of (44) is negative, we can writeso that π < k < 2π, which can also be expressed as

Now (29) takes the formwhere the sum is made over a set of cosine functions whose wave numbers are multiples of M, being possible to demonstrate (see the Appendix B) that the pattern formed by the density will have M peaks and M valleys. Note that, as then M ≥ 2. That is, the patterns will always have two or more peaks. Thus, those density distributions that have one peak will also be considered as SWPF.

In the case of B(t), Eq. (33) can be written aswhich establishes that, at least, the coefficients B 2(t) and B 3(t) need to be non-null, thus ensuring that B(t) ≠ 0. This means that the solution composed by the first three terms of the series in (47), not counting the coefficient B 0(t), represents the minimum configuration (MC) of non-null coefficients that can describe the formation of patterns.

Studying different configurations

Now let’s study the configurations obtained by considering a finite number of terms in (47) (not counting the coefficient B 0). We will analyze the behavior of the system given by the Eqs. (30) and (31) upon reaching a steady-state, obtaining the conditions under which patterns are generated.

One non-null coefficient

The density for this configuration will be given bywith which according to (31) implies that since B(t) ≠ 0. Then we can writebeing possible to determine the allowed values of M, B 0(t), and |B(t)| for a given combination of α and β, as shown in Table 1 . Note that (50) does not provide information about the sign of B(t), since the calculated quantity is .

Table 1

Values of M, B0(t), and |B(t)| for some combinations of α and β.

α=0.005			α=0.010
β=0.050			β=0.100
M	B0(te)	|BM(te)|	M	B0(te)	|BM(te)|
11	1.076	1.356	6	1.157	1.523
12	1.157	1.526	7	1.235	1.639
13	1.213	1.613	8	1.182	1.508
14	1.235	1.639	9	1.057	1.081
15	1.222	1.602
16	1.182	1.508
17	1.123	1.352
18	1.057	1.081
α=0.020			α=0.030
β=0.110			β=0.300
M	B0(te)	|BM(te)|	M	B0(te)	|BM(te)|
5	1.027	0.793	2	1.157	1.527
6	1.140	1.252	3	1.057	1.081
7	1.102	1.046

We must remember that the population density u(x, t) is a real and non-negative quantity, which will impose additional restrictions on the behavior of the coefficients. As B 0(t) is directly related to the size of the population, then its value is always non-negative. In the case of B(t), to ensure that u(x, t) ≥ 0, it is necessary that

When looking at Table 1, we note that the condition (51) is often unfulfilled. So the system composed of one non-null coefficient does not correctly describe the pattern formation. However, it is possible to obtain some results which will also be valid for more complex configurations. Thus, the relation (50) shows that to ensure that B 0(t) ≥ 0, it is necessary thatWith the increase in the number of non-null terms, the explicit form of B 0(t) would change, but this will not affect the validity of (52).

The expression for the coefficient B 0(t) calculated in (50) is practically equal to that given in (27), with the difference that now the value of k is fully defined. Therefore, we can use the following condition to determine the regions of the space (α, β) where pattern formation happens (B 0 > 1):Its validity will also be extended to other configurations since it does not depend on the coefficients. This relation allows the building of a phase diagram, as shown in Fig. 1 . For a given combination of α and β, we vary M between 2 and 50, keeping those values with which (53) is fulfilled. When at least one value fulfills the condition, we can say that the analyzed combination produces patterns. With this procedure, one obtains thatrepresents a necessary condition for pattern formation, thus confirming previous results [49], [50].

Fig. 1

Phase diagram obtained from (53).

From the results above, we can affirm that the study of systems simpler than the MC is equivalent to performing the perturbative analysis of the Eq. (18). With the increase in the number of non-null coefficients, a variation of B 0 should occur for each combination of α and β.

Two non-null coefficients

This time the density will bewith . At the same time ψ 2(t) ≠ 0, beingobtaining that

From (58) we get that and have the same sign. When the values of B 0(t) are limited to the following range:whereas if B 0(t) would have to comply that

The condition (61) is less restrictive, yielding values of B 0(t) higher than in the case of the system with one non-null coefficient. So we can affirm that, in the case of pattern formation, the sign of both and is negative, which leads to a positive sign of B 2(t).

For this configuration, it is also possible to build a phase diagram (not shown here) very similar to that of Fig. 1, but with higher B 0(t) values when compared with the previous case.

Three non-null coefficients

For the Minimum Configuration the density in the steady-state is given bywhere we removed the index t, since all operations will be done considering the steady-state.

We assume the fulfillment of (41), so thatUnlike previous configurations, now ψ(t) ≠ 0 for all the coefficients, so each of them will explicitly depend on the product of the other two.

Working with (32), (33), (36), and (37) we can write that Replacing now (66) into (65) yields Using the new expression for B 2 together with (66) in (64), leads to the following relation:which, after a few operations, is transformed into a quadratic equation for :where it was taken into account that ; ; with Finally, we can calculate as considering only the positive root of the discriminant.

Eq. (75) determines as a function of B 0 for a given combination of α, β, and M. With that it is already possible to calculate and . Knowing these values, we can use (67) for finding B 0 as a function of α, β, and M.

Before calculating the values of the coefficients, let’s make some observations about their signs. First, we can state, based on the results obtained when studying systems with a lower number of non-null coefficients, that and have a negative sign. Knowing this, and taking into account the hypothesis made on the behavior of the absolute value of the coefficients, we obtain, from (65), that the sign of B 2 is positive. It is now possible to establish, using the Eq. (66), that B and B 3 have the same sign, which again is undetermined since in the calculations we can only get .

As the absolute value of the coefficients falls rapidly with the increase of the index n × M, we can writewhich represents a sufficient condition for the MC to be able to correctly describe the pattern formed by some combination (α, β). That is, when (78) is fulfilled, the results provided by the MC should differ little from those obtained through a simulation that uses the same parameters, and that considers many more terms in the series given in (47).

Table 2 shows the values of the coefficients for some combinations of α and β. Note that B 2 is always positive, while the sign of both B and B 3 remains undetermined. The values of (B/B 3)2 are also presented, but in order to use the condition (78), we must first set a minimum value from which (B/B 3)2 ≫ 1. Here we will assume that if (B/B 3)2 > 30, then (78) is satisfied, and the MC will offer results similar to those of a simulation using the same values of α and β. This condition will be useful when comparing the results of the model with those provided by the simulations.

Table 2

Results of the Minimum Configuration for some combinations of α and β.

α=0.005						α=0.010
β=0.050						β=0.100
M	B0	|BM|	B2M	|B3M|	(BM/B3M)2	M	B0	|BM|	B2M	|B3M|	(BM/B3M)2
11	1.077	1.097	0.220	3.244·10^−3	1.144 · 105	6	1.166	1.640	0.542	4.938·10^−2	1.103 · 103
12	1.166	1.640	0.542	4.938·10^−2	1.103 · 103	7	1.316	2.037	0.922	2.190·10^−1	8.647 · 101
13	1.247	1.850	0.729	1.230·10^−1	2.262 · 102	8	1.687	4.440	3.214	1.620	7.514
14	1.316	2.037	0.922	2.190·10^−1	8.647 · 101	9	1.132	1.574	0.649	1.986·10^−1	6.279 · 101
15	1.368	2.223	1.137	3.500·10^−1	4.033 · 101
16	1.687	4.440	3.214	1.620	7.514
17	1.293	2.050	1.096	4.196·10^−1	2.388 · 101
18	1.132	1.574	0.649	1.986·10^−1	6.279 · 101
α=0.020						α=0.030
β=0.110						β=0.300
M	B0	|BM|	B2M	|B3M|	(BM/B3M)2	M	B0	|BM|	B2M	|B3M|	(BM/B3M)2
5	1.028	0.796	0.079	3.110·10^−4	6.546 · 106	2	1.166	1.640	0.542	4.938·10^−2	1.103 · 103
6	1.158	1.361	0.309	2.466·10^−2	3.048 · 103	3	1.132	1.574	0.649	1.986·10^−1	6.279 · 101
7	1.614	5.542	3.699	1.541	1.293 · 101

Degeneracy of states

Our model predicts the existence of several final states for the same combination of parameters. That is, each pair (α, β) presents some degeneracy, determined by the number of possible values of M. In general, degeneracy decreases with increasing parameter values, as shown in Table 3 .

Table 3

Degeneracy of states. For a fixed β, the number of admissible values of M increases as α decreases.

α∖β	0.050	0.100	0.150	0.200	0.250	0.300	0.350	0.400	0.450
0.090	-	-	-	-	-	-	1	-	-
0.080	-	-	-	-	-	-	1	1	-
0.070	-	-	-	-	-	1	1	1	-
0.060	-	-	-	-	1	1	1	1	1
0.050	-	-	-	1	1	1	1	1	1
0.040	-	-	-	2	1	2	1	1	1
0.030	-	-	2	2	1	2	1	1	1
0.020	-	3	3	2	1	2	2	2	2
0.010	6	4	4	4	2	4	3	3	3
0.005	8	7	6	6	4	5	4	5	4

We immediately start to ask which of the possible values of M will be the correct one, in the sense that a simulation could only present a fixed number of peaks in the density when reaching the steady-state. By having different final states, we could assume that the establishment of one of them is related to the initial conditions of the system. That is, changing the initial distribution of the population, it would be possible that the shape (number of peaks) of the final state also changes, for a given (α, β) combination.

If proven, state degeneracy would reveal a strong relationship between the initial conditions and the steady-state. But in the case of pattern formation, we expect that the final state is relatively independent of the initial conditions of the system [72], [73]. So degeneracy would be interpreted as a phenomenon that limits the combinations of α and β. That is, only those pairs that produce non-degenerate states (only one possible value of M) will be useful when trying to describe pattern formation.

High symmetry simulations

Besides the results corresponding to the analysis of the steady-state obtained previously, we can get additional information from the computational implementation of the system given by the Eqs. (30) and (31), comparing the model with the simulations.

The coefficient B 0(t) is described in (30) by a Riccati equation, while the others come from linear differential equations. Then we can use the method presented in Ref. [74] for numerically solve the system (30)–(31). In this way, we could explore the entire space (α, β), including those regions that are inaccessible through methods traditionally used to solve differential equations, such as the explicit Runge-Kutta methods.

So we are going to follow the time evolution of the first hundred terms of the series given in (29) (including B 0(t)) for different combinations of α and β. The initial values of the coefficients are calculated as choosing (unless otherwise specified)being whereas G is an arbitrary constant. The simulation continues until the steady-state is reached, which is assumed to be established when the conditionbecomes true. Here Δt represents the size of the time step used in the numerical integration.

Fig. 2 shows the time evolution of the coefficients B 0(t), B 1(t), and all those with an absolute steady-state value greater than 0.1. These values were obtained using (30) for B 0(t) and (31) for the others. The figure also shows the final density corresponding to each of the combinations (α, β).

Fig. 2

(a)–(d): Time evolution of the coefficients for different combinations of α and β. (e)–(h): Stationary values of the density. For each combination, the index M of the first non-null coefficient coincides with the number of peaks of the stationary density.

Note that for all the combinations shown. As B 0(t)L represents the total number of individuals, the value of B 0(t) must always be positive, which is confirmed by the simulations.

In Fig. 2(a) we see that the first non-null coefficient, excluding B 0, is B 9. All the others are null, except their multiples n × 9, and they become weaker and weaker, generating a figure with 9 peaks, as shown in Fig. 2(e). The same behavior is observed in (b), (c), and (d), changing only the number of peaks ((f), (g), and (h)). Altogether, we have the coefficients n × M, with M ≥ 2, which corresponds completely with the results obtained in Section 2.

Checking the degeneracy of the final states

As already stated, simulations will allow us to verify several results, such as the state degeneracy. To prove it, we will fix the values of α and β, varying only the density distribution at the beginning of the simulation.

Fig. 3 shows the results of the simulations made considering ; . We used initial distributions with 3, 4, 5, and 9 peaks, which resulted in final states with 6, 7, 8, and 9 peaks respectively, coinciding with the values of M reported in Table 2.

Fig. 3

Results of the simulations made considering ; and different initial distributions.

Degeneracy also appeared for the remaining pairs (α, β) reported in Table 2. Similarly, we studied several of the combinations presented in Table 3 as non-degenerate. For each of them we confirmed that the final state remained the same when the initial distribution was changed.

Thus, simulations confirmed the predictions made about the degeneracy of the states. There are pairs (α, β) for which the shape of the steady-state depends on the initial distribution of the system. At the same time, we found combinations whose final state is independent of the initial conditions. This last case seems to be more compatible with the principles of thermodynamics, since it reflects the fact that the system, during its evolution, eventually forgets its origins [73].

Comparing simulations with models

Knowing the combinations (α, β) that produce non-degenerate states, we can compare the results of the simulations with those provided by the models in Section 2. Fig. 4 shows the stationary density obtained at the end of the simulations, together with the one calculated from the coefficients provided by the MC. In most cases, the MC offers results very similar to those of the simulations.

Fig. 4

Stationary density for non-degenerate (α, β) combinations.

Meanwhile, Table 4 contains the values of the coefficients used to obtain the graphs of Fig. 4. When comparing, we noticed that precisely those coefficients that fail to fulfill the condition (B/B 3)2 > 30, generate patterns that do not match those provided by the simulations.

Table 4

Values of the coefficients provided by the MC in the case of non-degenerate (α, β) combinations.

α - β	M	B0	|BM|	B2M	|B3M|	(BM/B3M)2
0.030 - 0.250	3	1.575	4.325	3.065	1.277	1.147 · 1001
0.030 - 0.350	2	1.488	4.100	3.387	1.375	8.885
0.030 - 0.450	2	1.255	2.003	1.278	6.254·10^−01	1.025 · 1001
0.040 - 0.250	3	1.218	1.583	0.496	8.639·10^−02	3.358 · 1002
0.040 - 0.350	2	1.294	1.932	0.798	1.670·10^−01	1.339 · 1002
0.040 - 0.450	2	1.176	1.777	0.890	3.317·10^−01	2.868 · 1001
0.050 - 0.250	3	1.119	1.120	0.220	2.398·10^−02	2.181 · 1003
0.050 - 0.350	2	1.247	1.723	0.584	9.425·10^−02	3.343 · 1002
0.050 - 0.450	2	1.090	1.316	0.415	9.798·10^−02	1.803 · 1002
0.060 - 0.250	3	1.512	6.167	4.025	1.556	1.570 · 1001
0.060 - 0.350	2	1.196	1.496	0.405	4.976·10^−02	9.039 · 1002
0.060 - 0.450	2	1.022	0.655	0.083	8.173·10^−03	6.420 · 1003

Studying bacterial populations

The activities developed in Sections 2 and 3 have allowed us to make a broad assessment of the results obtained by solving the Eq. (18). Therefore, we are in an excellent position for trying to describe some experimental results.

We will highlight the work conducted by Anna L. Lin and collaborators, published in 2004 in the Biophysical Journal under the title “Localization and extinction of bacterial populations under inhomogeneous growth conditions”. This experiment (henceforward referred to as Ref. [41]) consists of a rectangular channel of length L, height h < L, and thickness E ≪ L, which is full of nutrients that can sustain a population of bacteria for a long time. The channel is illuminated with ultraviolet light, except for a region of length W covered by a protective plate (oasis), which moves with velocity v.

The species used in the studies was the E. coli RW 120. These bacteria die when exposed to intense ultraviolet radiation without undergoing mutations. With that, the population located under the oasis can grow normally. In the rest of the channel, the action of ultraviolet light inhibits the growth, eventually killing bacteria.

Each of the experimental runs started with the oasis located at the left edge of the channel. Then a small number of bacteria were inserted into the protected area. When (after several hours) the concentration of bacteria under the plate reached a saturation value, the oasis began to move at a constant speed v.

Between one execution and another, the parameter that usually changed was the speed of the oasis, whose values varied between cm/s and cm/s, so that the oasis could take between one and twenty days to reach the opposite edge of the channel. Other parameters were always the same, such as the channel dimensions (25 cm in length, 0.2 cm thick, and 2.5 cm height), the length of the oasis (3 cm, with thickness equal to that of the channel), and the intensity of the light source (10 Watt/m 2), consistent in a mercury lamp of 120 cm located 8.5 cm above the channel.

The vast majority of live bacteria remained close to the channel surface, while the dead ones took about two days to settle at the bottom. The channel had a measuring system attached that made it possible to obtain, every seven minutes, the value of the bacterial concentration near the surface throughout the analyzed space.

Using the Fisher equation to simulate the system

Besides the experimental activities, in Ref. [41] simulations of the system were performed using a generalization of the Fisher equation in the form [42], [46] which is basically the same as (7), with the difference that the growth term now has a coefficient that varies depending on the position of the oasis that moves with velocity v,being the position of the left edge of the oasis, while the value of ϵ is directly related to the intensity of the light.

The simulations performed in Ref. [41] managed to describe the general qualitative behavior observed for the bacterial density as a function of the speed of the oasis. That is, for relatively high velocities, bacteria cannot accompany the movement of the plate, so the entire population is exposed to radiation and eventually dies. At lower ones, part of the bacteria survives by moving together with the oasis.

Here the meaning of high or low speed is defined from the average speed value of the front of a bacterial population that grows, under favorable conditions, on a petri dish. Thus, oasis velocities close to are considered as high.

With the simulations, the experimentalists were able to make an extensive study of the parameters characterizing the population’s behavior, such as the velocity of the oasis, its length W, and the intensity of radiation. This assessment did not happen in the experiments due to their long duration.

Fig. 5 confirms several of the ideas exposed above. It presents the results of four simulations made in the present work at different velocities of the oasis, using (83) as a base, together with some experimental data reported in Ref. [41]. Concretely, we considered that cm/s, and cm/s. Other parameters, however, were set by us, such is the case of ϵ, that outside the oasis is considered to be zero, and the length of the oasis, fixed as being L exp the channel length, equal to 25 cm.

Fig. 5

Simulations performed using Eq. (83).

As in previous sections, here we used dimensionless constants and variables, which are obtained through the relations (12)–(17). In the case of the diffusion coefficient, the dimensionless value appears when dividing the experimental value by the quantity while the speed must be divided by aL exp. Thus we obtain that ; ; ; with .

In addition to the data shown in Fig. 5, we performed many other simulations using the Eq. (83). With them, it was possible to test several values of ϵ and W, always observing the same qualitative behavior described previously.

Although it does not appear explicitly, we can associate a value of the correlation length α to the Eq. (83). Working with (8) one obtains thatwhich yields cm, while .

In the calculations we assumed that [75] which allows us to rewrite (83) aswith . So that, at each point along the channel, the temporal evolution of the density u(x, t) follows a Riccati equation. For we had to solve a system of one hundred Riccati equations, for which we used the method presented in Ref. [74], considering .

Our simulations start with the left edge of the oasis located at the position being null its velocity. Initially, the population is allowed to grow only under the plate, which occurs employing a simulation similar to those described in Section 3, but considering . As (83) does not depend on β, at this stage we work with (18), using any value of β in the range 3.5α < β < W. In this way, we could test different distributions at the beginning of the oasis movement, but the results (not shown here) proved to be independent of them.

When the population density reaches a stationary value in the region protected by the oasis, it starts moving (from left to right) with speed v. During this stage we used homogeneous open boundary conditions, that is, if x < 0 or x > L. We specified the values of the density not at the boundary of the system, but beyond it. The justification for using this condition is that we are sure that the bacterial population does not exist outside the limits imposed by the channel. As the equations used in the simulations do not depend explicitly on ∂u/∂x, we don’t need to use Neumann boundary conditions. On the other hand, the situation considered in Sections 2 and 3 was highly symmetric, since the analyzed region was assumed to be one of the component cells of the world, being all the other cells the same, thus justifying the PBC. The experimental set is not as symmetrical, so that the use of PBC would induce unusual behaviors of the density, mainly at points near the boundaries.

When the right edge of the oasis reaches the position the simulation ends. So the time it takes for the oasis to reach the end of the channel is . Fig. 5 shows simulations made using the following oasis velocities: 0.25v; 0.5v; and 3.4v. For each of these velocities the density is shown in three different time instants: and .

A significant result obtained in Ref. [41] and confirmed here was that, although the parameters varied widely, the simulations failed to reproduce the double peak observed in the bacterial density for oasis velocities lower than . This discrepancy led the authors to suggest that the model used should be modified so that it takes into account the internal dynamics of various gene expression networks within individual bacteria. Thus, the observed oscillations could be the result of the variation in the levels of proteins present in bacteria as a response to environmental changes [76].

At this point, we can see a clear relationship between the experiment described and our model. On the one hand, the Fisher equation used in Ref. [41] is a particular case of Eq. (2), as shown in Section 2. While on the other hand, we already know that the solution of the Eq. (18) may generate different patterns. So we could assume that simulations based on our model should be able to describe the formation of several peaks in the bacterial density for certain combinations of α, β, v, and W.

Considering non-locality

Now we will use (2) to perform simulations. Let’s assume that being possible to rewrite (2) aswhich differs from (83) due to the inclusion of the non-local kernel dependent on β, so that this parameter will be responsible for any difference observed in the results.

Considering (86), the kernel given in (9), and using the trapezoidal rule to solve the integral, we obtain, after conveniently grouping the terms:where dimensionless constants and variables were also inserted. Here whereas that N is the integer closest to 2β/Δx. Again we have a Riccati equation, so the numerical solution will be obtained by applying the method presented in Ref. [74], with the same conditions and values of the parameters used previously.

The information shown in Fig. 6 is similar to Fig. 5, but in the case of Fig. 6 the simulations used the Eq. (91), with . We immediately observe the presence of double peaks in density for velocities lower than v, what didn’t happen in Fig. 5. The best definition of the double peak takes place for oasis speeds close to 0.25v, which is in total correspondence with the experimental values reported in Ref. [41], where the appearance of more than one crest occurred for v/v ≈ 0.23. As the velocity of the oasis increases, the intensity of one of the peaks decreases, remaining only one when v ≈ v.

Fig. 6

Results of the simulations carried out considering the parameter β.

In Fig. 7 , it is possible to appreciate in more detail the temporal evolution of the system when the speed of the oasis is . Initially, the population in the region protected by the oasis increases so that its distribution is asymmetrical, with the peak shifted towards the left edge of the oasis. The amplitude of this peak decreases as another one appears at its right. When the peaks are almost equal. From that instant, the crest on the right will have a larger size. The waveform that moves with the oasis will preserve these overall features until the end of the simulation.

Fig. 7

Temporal evolution of the system with oasis speed and .

The validity of some results obtained in Sections 2 and 3 was maintained, even when the equation and boundary conditions were not the same. Such is the case of the value β from which the double peak begins to appear. When the double peak appears for β > 0.17, which is in correspondence with Table 3. The aforementioned confirms the generality of our model, placing β as the most important parameter for pattern generation.

Finally, we can say that to observe more than one peak when making simulations using the Eq. (90), besides the condition β > 3.5α, it is also necessary thatsince the region where the growing conditions are favorable (similar to those used in Section 3) should be large enough. Condition (92) was proven in the simulations, thus justifying the value of the length of the oasis in Figs. 6 and 7, which is larger than the one used in Ref. [41] ().

Understanding the parameter β

The analysis of Fig. 5, Fig. 6, Fig. 7 reveals that the parameter β is responsible for the appearance of a double oscillating peak in bacterial density during the movement of the oasis through the channel. But there is still no expression that relates β to experimental quantities.

Based on Eq. (85), we can think of the quantity as being the mean square displacement performed in the time unit by any of the bacteria, under favorable growth conditions. This definition has a statistical character due to the dependence on the diffusion coefficient). It establishes the connection between one of the parameters of our equation and an intrinsic property of the system; namely, the displacements of bacteria when the growth rate is the same across space.

We could assume that β relates in a similar way to the studied system, but this time the displacements would be associated with non-homogeneous growth conditions. That is, all those factors that provide unequal growth will combine so that, besides the normal one, there is an additional movement component in response to adverse external conditions, whose origin lies in internal processes that take place in bacteria.

Searching in the literature [76], [77], [78], we found that many organisms use the so-called circadian clocks as sensors of environmental activity, which allows them to react appropriately to external stimuli. Most of these clocks use intracellular gene expression networks composed basically of positive and negative regulatory elements whose interaction gives rise to oscillations of gene expression. The positive component activates genes associated with the clock and promotes the manifestation of the negative factor, whose objective is to decrease the concentration of the positive one. The cycle closes with the degradation of the negative element and the reappearance of the positive agent.

A crucial feature of the clocks described above is that they can maintain a constant period [79], even under the influence of external and internal fluctuations. Among them, we have temperature variations, chemical reactions that happen randomly within bacteria, abnormal concentrations of some chemical species. Thus, a significant group of processes can occur at the correct time.

But there must be some limit on the amplitude of the fluctuations below which the clocks can operate normally. So one might ask, what happens beyond that limit? At this point, we can return to the experiment carried out in Ref. [41], where part of the space was under the influence of very intense radiation. This radiation represents a modification of the environment so significant that it irreversibly affects the correct functioning of the clocks within the bacteria, eventually causing their death.

It is also logical to assume that before the collapse of the bacterium occurs, some mechanism is activated that warns of the impossibility of maintaining the correct performance of vital functions under the current environmental conditions. That is, an alert message is issued asking the bacterium to run for its life. Thus, the appearance of an additional component in the movement of bacteria would represent, in the case analyzed, a desperate response, aiming to find better living conditions. Note that, by its nature, the new component will always be larger than the normal one, which is in correspondence with our model, where patterns arise, in general, when β > 3.5α.

Then we could propose, by analogy with (85), the following relationship for the parameter β:being D the emergency diffusion developed by bacteria subjected to unfavorable conditions. Here D represents an average of the squared displacements of bacteria. But the functional dependence of each of these displacements is, in general, rather complicated, since it will be related to internal and external factors acting on each bacterium.

Among the internal factors, we find the physical-chemical state of the bacterium, determined by the correct functioning of the internal clocks and by the amount of accumulated energy (obtained from the consumption of nutrients) at the moment when it decides to escape. Thus, a bacterium whose systems are functioning correctly when the environment becomes hostile must travel greater distances than another bacterium of the same species previously damaged. Similarly, the higher the accumulated energy, the higher the kinetic energy at the time of escape, which will result in a greater distance covered in the unit of time.

As an external factor, we have the spatial location, since a bacterium positioned in an area with relatively low population density will collide less with its neighbors, being able to travel greater distances than another bacterium located in a high-density area. However, at low concentrations, the bacterium could choose the wrong path and end up dying, since communication mechanisms such as Quorum Sensing [80], [81], [82], [83], [84], [85], [86] are not sufficiently stimulated.

The elements listed above were analyzed individually, but in a real situation, they should interfere, giving way to a complicated movement that would only be accessible through the performance of specially prepared simulations.

Finally, it is possible to say that the parameter β represents a macroscopic measure of changes that occur at an individual level in the concentrations of proteins in bacteria, as a product of significant variations in the environment. Then the magnitude would represent the mean square displacement performed in the unit of time by bacteria under unfavorable growth conditions.

Going back to degeneracy

Now we will refer to a result obtained in Sections 2 and 3, which is, to some extent, intriguing: the existence of degenerate states for some (α, β) combinations. Due to the marked dependence on the initial conditions, these states were discarded as viable candidates to describe the patterns formed by the population in response to sudden changes in the environment. But the question of whether it would be possible to associate them with other systems or processes remained unanswered.

According to Table 3, for a fixed β, the degeneracy increases with the decrease in the value of α. This parameter is closely related to the diffusion coefficient, as shown in (85). A reduction in the value of the diffusion coefficient leads to a rise in the time required for the system to reach the steady-state (relaxation time). So that there could be pairs (α, β), mostly degenerate, whose relaxation times are much longer than those required by the gene expression networks to adapt to fluctuations in the environment.

But the time scale of these pairs could be compatible with the times associated with processes such as mutation, whose manifestation in the characteristics of the population is, in general, more delayed. This suspicion is reinforced by Fig. 8 , which shows the results of simulations carried out following the specifications of Section 3, corresponding to the combination (α, β) = (0.01, 0.10). We used two different initial conditions, so next to each of them, we also show the final distribution of the density, as well as the time evolution of the coefficient B 0, whose value is directly proportional to the size of the population.

Fig. 8

Influence of the initial conditions.

The coefficient B 0 in Fig. 8(f) presents during its temporal evolution what we could call a metastable state, during which its value remains practically constant, subsequently occurring the transition to the final stage. At this point, we could associate the initial conditions with different mutagenic agents affecting the population.

The metastable state could represent the stage during which the entire population is affected by the agent, thus ensuring the presence of the mutation in the following generations. This mutation would manifest itself through a collective response that varies according to the agent (dependence between initial and final states).

However, the previous analysis is still at the level of speculation. It serves to establish an interest in further study of the residual information provided by our model, assuming that it can describe other real processes.

Conclusions

The work done so far can be divided into two large blocks. On the one hand, the use of RDE to determine the conditions under which patterns appear. On the other hand, activities aimed at both describing experimental results and establishing the physical meaning of parameters associated with the model developed in the first block.

Initially, we presented a set of equations based on non-local operators for studying the spatial patterns formed by a population. The solution took the form of a Fourier series, which allowed an optimal use of the advantages offered by the periodic boundary conditions. Thus, starting from a handful of assumptions on the behavior of the expansion coefficients, we could establish the main features of the patterns in the steady-state, which can be summarized as follows:

* The parameters α and β determine the existence of a phase space, characterized by regions with and without patterns;

* The patterns formed will have M peaks and M valleys, being M the index of the first non-null coefficient of the series given in (29), with M ≥ 2;

* There are (α, β) combinations, which lead to the formation of degenerate states, where the system seems not to forget its initial conditions.

These results were confirmed by the simulations in Section 3, allowing to determine those (α, β) combinations that could be useful to describe patterns.

The second block was already partly built, since previous works established a relationship between the parameter α and the diffusion coefficient [49], [50], [87]. Here we extended the idea to the case of β, proposing the relation (93). Now the movements would not be random, but they would have the objective of making the bacterium manages to leave that region of space that did not offer conditions for the maintenance of life. In other words, we obtained that, besides the random component, bacteria would develop a second one in their movement, more organized, as a response to unfavorable environmental conditions.

The physical interpretation of the parameter β is indissolubly related to the success we had in adapting our model to the conditions of the experiment. In that process, we lost the high symmetry of Sections 2 and 3, despite which several results remained valid. Finally, the simulations made in Section 4 showed that our model, unlike others previously used, can reproduce some experimental results, such is the case of the double peak in the bacterial density during the movement of the oasis through the channel.

This last result is very significant, as it demonstrates the usefulness of our model and leaves the door open for future research. And several activities could be developed. Firstly, we have that the present work focused on the study of a population in the one-dimensional case, which offered quite good results. So that a natural extension would be the study of the model in two dimensions, considering one or several species interacting in the same space, including in this way prey-predator interactions.

When introducing the RDE, it became clear that the form of the source term would be of great importance, since it is directly related to the characteristics of the system. In the case of our model, it is the form of non-local kernels that has a significant weight. We chose the one given in (9) to facilitate mathematical work. So it would be interesting to use another functional form (for example, something like an exponential decay) and observe the changes that can happen.

The simulations of Section 4 considered nutrients as equally distributed throughout the space, being intense radiation, the only stimulus affecting the behavior of bacteria. Nevertheless, it is highly probable that there is a heterogeneous spatial distribution of nutrients, giving rise to phenomena such as chemotaxis (movement of organisms towards or away from a chemical) [88], that can make the diffusion coefficient is dependent on the position and time [4]. These variations in the diffusion coefficient could cause additional oscillations in density, which would be visible at shorter oasis lengths than those predicted by our model, which is what happens in the experiment. So incorporating chemotaxis could produce more realistic results.

CRediT authorship contribution statement

Orestes Tumbarell Aranda: Conceptualization, Methodology, Software, Validation, Formal analysis, Data curation, Writing - original draft, Writing - review & editing, Investigation. André L.A. Penna: Conceptualization, Methodology, Writing - original draft, Writing - review & editing. Fernando A. Oliveira: Conceptualization, Methodology, Writing - original draft, Writing - review & editing, Supervision.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.