Turing Instability in an Economic-Demographic Dynamical System Can Lead to Pattern Formation on Geographical Scale

Spatial distribution of the human population is distinctly heterogeneous, e.g. showing significant difference in the population density between urban and rural areas. In the historical perspective, i.e. on the timescale of centuries, the emergence of the densely populated areas at their present locations is widely believed to be linked to more favourable environmental and climatic conditions. In this paper, we challenge this point of view. We first identify a few areas at different parts of the world where the environmental conditions (quantified by the temperature, precipitation and elevation) are approximately uniform over thousands of miles. We then examine the population distribution across those areas to show that, in spite of the homogeneity of the environment, it exhibits a clear nearly-periodic spatial pattern. Based on this apparent disagreement, we hypothesize that there exists an inherent mechanism that can lead to pattern formation even in a uniform environment. We consider a mathematical model of the coupled demographic-economic dynamics and show that its spatially uniform, locally stable steady state can give rise to a periodic spatial pattern due to the Turing instability. Using computer simulations, we show that, interestingly, the emergence of the Turing patterns eventually leads to the system collapse.


Introduction
Fast growth of the global human population has long been regarded as a major challenge that faces the mankind [35,22,52,53]. Presently, this challenge is becoming even more serious than before, in particular because many natural resources are estimated to deplete before the end of this century. The increasing population pressure on the agriculture and on the ecosystems and the environment more generally is predicted to result in worldwide food and water shortages, pollution, lack of housing, poverty and social tension. The situation is exacerbated by the global climate change as considerable areas of lend are predicted to be flooded and hence taken out of humans use. It is widely believed that, unless alternative scenarios of sustainable population growth and social development are identified and implemented, the mankind is likely to experience stagnation or even decline [36].
Population growth in time is complemented with the population dynamics in space. Population distribution over space is hugely heterogeneous for a variety of reasons, to mention the climate, the history, and the economy just as a few. The spatial heterogeneity may result in significant migration flows that in turn can have a significant feedback on the local demography and the population growth. On a smaller scale of individual countries and states, understanding of the factors affecting the population distribution in space is needed to ensure adequate development of infrastructure, transport, and energy network. Poorly informed decisions are likely to result in overcrowding and social problems in urban areas and/or lower quality of life in rural neighbourhoods.
Identification of scenarios of sustainable population growth and social development on various spatial and temporal scales requires good understanding of the relevant processes and mechanisms that affect both the population growth and the population distribution. Arguably, such understanding is unlikely to be achieved without a well-developed theory and the corresponding mathematical/modelling framework. Indeed, mathematical models of human population dynamics (e.g. see [13,31]) has a long history dating back to the 17th century [29]. Over the last few decades, the need for an adequate and efficient mathematical theory of the human population dynamics has been reflected by a steady growth in the number of studies where problems of demography were considered using mathematical models, tools and techniques; see Fig. 1.
In this paper, we use mathematical modelling to address the phenomenon of heterogeneous spatial population distribution. Heterogeneity of geographical features (mountains, forests, rivers, etc.) and natural resources (e.g. coal, iron and copper ore) are commonly accepted as factors leading to the demographic and economic heterogeneity. However, the question that we ask here is -is this natural heterogeneity the only underlying cause, or can there be another and perhaps more general principle responsible for emergence of heterogeneous population distribution?
In order to answer this question, we first revisit available data on the population density over a few areas in different parts of the world to show that, in all cases, the population distribution exhibits a clear nearly-periodic spatial pattern in spite of the fact that the environmental conditions are relatively uniform. We then consider a novel model of coupled economic-demographic dynamics in space and time and endeavour to use it to simulate the spatial population distribution. The model consists of two coupled partial-differential equations of reaction-diffusion type. We show that the emergence of spatial patterns appears to be possible as a result of Turing instability. By relating the model predictions to the data on the human population density, we argue that the heterogeneous population distribution observed across different countries in different continents may have been caused by endogenous rather than exogenous factors, i.e. may have appeared due to intrinsic Turing instability of the corresponding economic-demographic dynamical system.

Real-world examples
In many countries, the population distribution over space is distinctly heterogeneous, e.g. urbanized areas with a high population density alternate with rural areas with a low population density. Apparently, spatial variation in geographical and climatic factors can play a significant role in shaping the population distribution. Since our main hypothesis in this paper is the existence of a dynamical mechanism that may lead to the formation of heterogeneous population distribution regardless of the geographical heterogeneity, in our search for the real-world examples we focus on the cases where the environment may be regarded as relatively uniform. The environmental properties that we consider here as proxies for the environmental heterogeneity are the elevation, the annual mean temperature, and the annual mean precipitation. A brief overview of the several relevant cases is given below.

Canadian Southern Region
Canada is a scarcely populated country and the majority of Canadian population live in the narrow band (approx. 160 km) along the USA border, see Fig. 2a. The distribution of the environmental properties across the country is highly heterogeneous, in particular in the South-North direction, ranging from temperate climate in the South to the rather extreme polar climate in the North. However, the magnitude of climatic variation in the East-West direction is much smaller (see Fig. 2c), at least over the span between the Atlantic coast and the Rocky Mountains where the annual mean temperature varies just within 2-3 • C (contrary to about 20 • C in the South-North  [45]. One can see that, to the east of the Rocky Mountains, the elevation along the Canada-USA border is approximately uniform. (c) Annual mean temperature map of Canada in 2016 [12]. It is readily observed that the temperature does not vary much along the southern border. (d) Precipitation map of Canada [25]. The amount of precipitation does not vary much along the border, except for the extreme West. direction). A similar observation applies to the elevation and the annual mean precipitation; see Figs. 2b and 2d, respectively. We now focus our analysis on the curved narrow corridor along the border (see the green line in Fig. 2a) where the environmental conditions are relatively uniform but the population distribution is not. Figure 3a demonstrates how the population density inside the band varies in space along the East-West direction. Interestingly, the distribution exhibits three maxima with approximately equal spacing of 700 miles. We therefore regard it as a periodic spatial distribution. Note that this pattern is persistent over time: a similar periodic-like structure is observed for different years (not shown here for the sake of brevity) starting from at least late 19th century.
In order to reveal how strong is the effect of environmental properties on the population distribution, we now perform the pairwise correlation analysis between the population density and

South-Eastern Australia
As another relevant example of the heterogeneous population distribution in an approximately uniform environment, we consider South-Eastern Australia. As well as Canada, Australia is a  as precipitation and temperature are approximately uniform (see Figs. 4b-d), e.g. the variation in the annual mean temperature is just a few degrees (compared to more than 30 • C over the continent as a whole). The stripe includes the Great Dividing Range and the Australian Alps, which therefore accounts for a significant variation in the elevation. In spite of the relatively uniform environment (apart from the elevation, its effect is discussed below), the population distribution along the strip is clearly heterogeneous with the population density varying several times between the more dense areas and the less dense ones; see Fig. 5a. Interestingly, it exhibits a nearly-periodic pattern where the three maxima are approximately . In all cases shown in panels (b-d), the best-fitting straight line is drawn by maximizing R 2 ; for details, see Table 1. equally spaced by about 700-850 kms. An immediate intuitive explanation of the heterogeneous population distribution can be sought in the heterogeneity of the environmental properties. Correspondingly, we look into the effect of the environmental factors more carefully by considering the correlation between each of the three factors chosen above and the population density. Figures 5b and 5c show the scatterplots of the population density in Australia vs the mean annual temperature and the mean annual precipitation, respectively. In both cases, the straight line shows the best-fitting of the data to maximize R 2 ; the corresponding values of R 2 are shown in Table 1. It is readily seen that in both cases R 2 is quite small. We therefore conclude that the climatic variation is unlikely to be the factor that defines the spatial distribution of the population. Now we recall that the study area includes the mountain ranges and exhibits considerable variation in the elevation. The question hence arises as to whether that can be a relevant factor. However, we first notice that the vast majority of the Australian population leaves at the elevation below 250 meters; see Fig. 5d. We then perform the correlation analysis by looking for the best-fitting straight line in the scatterplot of the population density vs the elevation. The corresponding value of R 2 (see Table 1) appears to be very small. We therefore rule out the elevation as a factor affecting the heterogeneous spatial population distribution along the South-East coast of Australia.

Mongolian Grassland
Mongolia, a Central-Asian country situated between China in the South and Russia in the North, has an elongated territory that extends from East to West for about 2400 km. It is the most sparsely populated country in the world. South of Mongolia is occupied by the Gobi Desert, which is barely populated at all due to the harsh climate and lack of resources. The majority of Mongolian three million population live in grasslands, which is located in the North of the country. In order to reveal the features of the spatial population distribution as is needed in the context of this study, we focus on the densely populated narrow corridor located along the latitude at 47.7 degrees North; see the black line in Fig. 6a. Interestingly, we readily observe that, as well as in the two previous cases, the population distribution in the East-West direction exhibits a periodic-like pattern (Fig. 7a). The three distinct peaks are separated by 700 and 900 kms intervals.
Variation of the environmental properties (cf. Figs. 6b-d) along the latitude is considerably less than in the North-South direction. However, it appears to be larger than it is in the cases of Canada and Australia, e.g. the annual mean temperature varies over about 10 • C and the annual mean precipitation from 50 to 350 mm/(m 2 ·year). Also the elevation varies over about 1500 meters, which is somewhat less than in Australia but larger than in Canada (where our analysis did not include the Rocky Mountains).
In order to reveal whether the variation of the environmental properties has any significant effect on the distribution of the population, we now perform the pairwise correlation analysis. The scatterplots of the population density vs the mean annual temperature, mean precipitation and the elevation are shown in Fig. 7b, c and d, respectively. The straight line is the best-fitting linear function; the corresponding values of R 2 are given in Table 1. Apparently, the correlation between the population distribution and the environmental factors is very weak. We therefore conclude that the nearly-periodic pattern clearly seen in the population distribution is unlikely to be caused by the environmental conditions.

Multiple linear regression
In the above, we have shown that the spatial distribution of the population is unlikely to be affected, not to any considerable extent, by any single environmental property such as the mean annual temperature, mean precipitation or the elevation. However, generally speaking this does not rule out a possibility that a certain combination of those three factors may have a much stronger effect. In order to check this possibility, we applied the multiple regression: where y is the population density, x 1 , x 2 and x 3 are, respectively, the average annual temperature, the annual precipitation and the elevation. Model (1) was applied separately to the data for each of the three countries. The results are shown in Table 2. We readily observe that the joint consideration of the three environmental factors does not real to a stronger correlation. Small values . In all cases shown in panels (b-d), the best-fitting straight line is drawn by maximizing R 2 ; for details, see Table 1.
of R 2 indicate that variation of population density is only to a small amount explained by the geographical and climatic properties. Thus, we have examined three areas in three different countries chosen from three different continents to reveal that, in all three cases, the population distribution over an area with relatively uniform environmental conditions exhibits a clear spatial periodicity. Having considered the correlation between the population density and the main environmental properties, we have shown that the correlation is very weak and hence the nearly-periodic pattern is unlikely caused by the effect of the environmental factors. Note that the three considered countries are vastly different in term of their average climate, history, ethnicity and culture. This leads us to assume that there can be a generic mechanism resulting in the emergence of the observed spatial pattern. We further assume that this is a dynamical mechanism originated in the nonlinear interaction between the human demography and the distribution of resources or wealth. The corresponding mathematical model is considered in the next section.

Mathematical model
In order to describe the dynamics of the human population, we use the simple, "conceptual" economic-demographic model earlier developed in [55]. The model quantifies the state of the human society at a given location in space x at a given time t by two state variables, the population density p(x, t) and the concentration of wealth u(x, t). Note that, whilst due to its meaning p ≥ 0, variable u must not necessarily be non-negative; negative values of wealth can be regarded as debt.
In the baseline 1D case (which is relevant in case of the population distribution in a narrow stripe, cf. the examples in the previous section), the model consists of two partial differential equations of reaction-diffusion type: (where we neglect, for the sake of simplicity, possible effects of cross-diffusion [55]). Here the first term in the right-hand side of Eq. (3) accounts for the population movement in space, that we assume can be considered, at least over certain spatial and temporal scales, as random [14,46,58] (for a detailed discussion of the "bugbear of randomness" see [46]) and can be described mathematically as standard Fickian diffusion. The diffusion term in Eq. (2) describes local wealth redistribution due to the economic activities such as trade and investments, and/or taxes.
In order to specify the population growth, we consider the reaction term in Eq. (3) in the following form: where the first term in the right-hand side describes the reproduction rate of the population. We therefore consider it to be the standard logistic population growth with the fertility rate α and the carrying capacity K. The second term describes mortality, σ being the mortality rate. Both the carrying capacity and the mortality rates depend on the wealth. For the mortality rate, we consider it to be a monotonously decreasing function of wealth. It takes into account the general observation that, on average, the mortality rate is lower for rich people, e.g. due to access to better health services and/or healthier life style [19]. In particular, there is evidence that in the USA wealthier people tend to live longer [10]. More specifically, we consider the following generic Monod-type parametrization: where c 0 , σ 0 and σ 1 are positive parameters, σ 0 > σ 1 . In order to parameterize K(u), we first recall that it describes an equilibrium population density allowed by the availability of resources [42]. When the resource -in our case, wealthbecomes scarce, the carrying capacity goes to zero, K(0) = 0. Therefore, for small u, K(u) is an increasing function. However, when the resource (wealth) becomes a plenty, K(u) seizes to be monotonic. There is a certain cultural shift between the low-income and high-income society groups [43]. For people with a low income wealth is the main limiting factor, whilst for people with a high income it is not necessarily so. In particular, low-income people typically live in urban areas, and hence at relatively high population density, e.g. due to their dependence on public transport and other public services [27]. With an increase in income, people tend to move to less densely populated areas such as a rich suburb or a large private estate. Correspondingly, we consider the carrying capacity K(u) to be an increasing function of wealth for small u but decreasing function for large u, tending to a small value (ultimately, to zero) as u tends to infinity. More specifically, we consider the carrying capacity in the following form: where a 2 and c 2 are positive parameters. In order to specify the reaction term in Eq. (2), we first write it as follows: where W and S are the rates of the wealth production and consumption, respectively. Production of wealth is often described by the Cobb-Douglas production function which in the simplest case can be written as [30,44] where L is the labour, Q is the capital and M is the available natural resource, a positive coefficient b is a measure of technology, ν, β and γ are positive constants [30,44]. We assume that the natural resource is not a limiting factor so that M can be kept as constant. We further assume that capital Q is a function of wealth, Q = f (u), and labour is a function of the population density, L = g(p). Equation (8) then takes the following form: Due to their meaning, it is reasonable to assume that f (u) and g(p) are increasing functions with saturation. Correspondingly, we choose them in the generic form as the Monod function: where a 1 , c 1 and c 2 are positive parameters. For the wealth consumption S, we assume it to be the result of two processes, i.e. due to the depreciation (in particular in case of buildings, machinery, etc.) and due to the consumption of the goods and products by the people. For the depreciation, we assume it to be a linear process with a constant rate a. The rate of the individual (per capita) consumption, say c, can be described by the Keynes linear consumption function, c = r + sy, where y is the per capita income and r and s are positive coefficients. Assuming additionally that average income is proportional to the wealth, we arrive at the following expression: S(u, p) = au + (r + su)p.
From (4-11), we thus obtain the following expressions for the reaction terms:

A glance at the nonspatial system
We begin with a brief look at the properties of the nonspatial counterpart of the reaction-diffusion system (2-3), which is given by the following equations: where functions G and F are given by Eqs. (12)(13). System (14) was studies in detail in [59].
Here we only briefly revisit some of its properties, to the extent that is needed for the goals of this paper. The phase plane of system (14) is shown in Fig. 8. It is readily seen that the origin is a steady state; a closer look reveals that it is stable node. Inside the first quarter of the phase plane, i.e. for u > 0 and p > 0, the F -isocline is a convex closed curve (loop) and the G-isocline is an upwardconvex, dome-shaped curve. Depending on the relative position of the isoclines (and hence on the parameter values, see [59] for details), the number of the positive steady states can be anywhere from 0 to 4. Therefore, in a general case system (14) can exhibit a rich, multi-stable dynamics and a complicated bifurcation structure where positive states can emerge or disappear. A typical case corresponding to four positive steady states is shown in Fig. 8a.
A case where the relative position of the isoclines allows for only two positive steady states is shown in Fig. 8b. For these parameters, A is a saddle point and B is a stable focus. Interestingly, a closer look reveals that even in this case the phase plane has a complicated structure; see Fig. 9. There are two attractors: the stable node (0.0) and the stable focus B, so that the system is bistable. The attraction basin of stable focus B is bounded by an unstable limit cycle (shown by brown colour). Trajectories that start close to the limit cycle from inside will in the large time limit approach the stable focus; an example is shown by the red curve. We menton here that the eigenvalues of the system linearized in the vicinity of stable focus B have very small real part (for the parameters of Fig. 9, λ 1,2 = −0.0000895 ± 0.460312i) so that the trajectory approaches the steady state at a very low rate. Trajectories that start outside of the limit cycle eventually approach the stable node (0, 0) except for the special trajectory (the blue curve) that is a part of the stable manifold of saddle point A; an example is shown in Fig. 10.

Turing instability conditions
We now consider the properties of the spatially-explicit system (2-3) with the reaction terms given by (12)(13): Equations (15)(16) are complemented with the Neumann 'zero-flux' boundary conditions: Since our study is motivated by the existence of periodic spatial patterns, see Section 2, we are particularly interested in the possibility of Turing instability and the corresponding pattern formation. Turing instability is the property of a nonlinear reaction-diffusion system where a steady state that is stable in the corresponding non-spatial system can, under certain parameter constraints, become unstable in the spatial system with respect to periodic perturbation with a certain wavelength [51]. The brown curve shows the unstable limit cycle. Note that the green curve leaves the vicinity of steady states A and B to eventually go to the stable node (0, 0) which is the attractor for the rest of the phase plane (except for the part of the plane inside the limit cycle). The arrows show the direction of the phase flow. The red curve shows a trajectory that starts close to the limit cycle and eventually converges to the stable focus.
Let (ū,p) is a steady state of the nonspatial system and J is the Jacobian matrix evaluated at this steady state: where the subscript denotes the corresponding partial derivative, for instance F u = ∂F (u, p)/∂u. We require that the steady state is stable, so that the following conditions hold: (a) tr(J) < 0, e.g. see [21]. In the spatial system, for stability of the corresponding uniform steady state p(x, t) ≡p and u(x, t) ≡ū with respect to a periodic perturbation with the wavenumber k, conditions (19) change to (a) tr(J k ) < 0, where The instability occurs if one of the conditions (20) is broken. It is readily seen that condition (20a) holds for any k. Therefore, the instability can only occur if there is a range of values of k that satisfy the following inequality [21,34]: where the characteristic function F (k 2 ) = det(J k ). Taking into account (21), F (z) appears to be a square polynomial, so that inequality (22) is equivalent to [21,51]: where parameter D = D p /D u is the ratio of the diffusion coefficients. In its turn, it appears that a necessary condition for (23) is that F u and G p must be of different sign. Consider F u > 0 and G p < 0; in this case, u is called the "activator" and p the "inhibitor" [21]. Then another necessary condition for (23) is D > D cr > 1 where D cr is a certain critical value that depends on the parameters in the reaction terms [21,51]. Now we consider how the generic relation (23) between the system's feedbacks works in the case of our model (15)(16). Given the complexity of the bifurcation structure of nonspatial system (13)(14), see Section 4 and [59] for more details, a comprehensive study addressing the Turing instability of all (stable) steady states over the entire parameter range does not seem possible. We therefore concentrate on the specific yet instructive case where there are two positive steady states, a saddle and a stable focus (cf. Figs. 8b and 9), in particular to investigate whether the Turing instability may occur for stable steady state B = (ū,p).
As a starting point, we consider the following set of parameter values: a 1 = 170, c 1 = 5, c 0 = c 2 = 1, a = 0.01, r = 7.5, s = 7.5, a 2 = 20, c 3 = 0.5, σ 1 = 0.05, σ 0 = 190, α = 10. The corresponding steady state values areū = 0.69 andp = 0.00087 and the Jacobi matrix is Therefore, at this steady state wealth acts as the activator and population as the inhibitor: the relation similar to the classical resource-consumer system. For these parameters, the critical value of the diffusivity ratio is readily obtained as D cr = 59746. Figure 11 shows the function F (z) for a subcritical case < cr where the steady state is stable (as F (z) > 0 for any z and condition (20b) holds for any k) and a supercritical cases > cr where the steady state is unstable with respect to perturbation with the wavelength from the interval where F (z) < 0 and hence condition (20b) is broken.
For the above parameter set, the critical ratio of the diffusion coefficients is very large, which may rise doubts whether it is at all realistic in terms of the real-world dynamics. Therefore, now we are going to consider how the critical relation responds to changes in the parameter values and whether it can be diminished. Indeed, it appears that D cr is rather sensitive with respect to the variation of some of the model parameters; examples are shown in Fig. 12. We have found that by varying α, a 1 and s, the critical ratio can be made as small as D cr = 98.5 (obtained for parameter values c 1 = 5, c 0 = c 2 = 1, a = 0.01, r = 7.5, s = 16, a 2 = 20, c 3 = 0.5, σ 1 = 0.05, σ 0 =  Fig. 11. 190, α = 9.63, a 1 = 168, the corresponding steady state values areū = 0.59 andp = 0.001).
We mention here that a further reduction of D cr does not appear to be possible: for instance, a further decrease in α (as in Fig. 12a) or a further increase in s (as in Fig. 12b) make the steady state unstable.

Spatiotemporal dynamics: numerical results
In this section, we consider the spatiotemporal dynamics of system (15)(16) that arises as a result of the Turing instability. Note that the fact that the steady state becomes, in a certain parameter range, unstable with regard to spatially heterogeneous perturbations is established analytically (see Section 5) and hence, as such, do not require any confirmation (e.g. by simulations). However, the mathematical analysis of the instability is based on the system linearization in the vicinity of the steady state and thus is limited to the time interval when the deviation of the solution from the steady state is small. That rises a question what can be the dynamics at the later time, after the deviation from the steady state becomes large enough to be affected by the nonlinearity of the system. Turing instability is known to often lead to the formation of a stationary spatially-periodical pattern [37], however more complicated dynamics can occur too [18,33]. In order to make an insight into the above question, the reaction-diffusion system (15)(16) is solved numerically by finite-differences using the following initial conditions: with the size of the spatial domain L = 120. The diffusion coefficients are chosen as D u = 1 and D p = 100, and the values of the reaction parameters are given at the end of the previous section.
At the boundaries of the domain, the zero-flux conditions (17) are used. For these parameter values, D p /D u > D cr = 98.5 so that we expect that a small initial perturbation of the steady state leads to the emergence of a spatial pattern. This is indeed seen in the numerical simulations; see Fig. 13. At an early stage of the dynamics, the initial conditions (25) fast evolve (over t ∼ 100, not shown here) to a nearly stationary periodic spatial pattern, which then remains almost unchanged over a considerable time: until t ≈ 17200, see Fig. 13a). The spatial distribution then starts evolving at a much faster rate, first showing a considerable (nearly three-fold) increase in the spatial variability, see Fig. 13b, and then developing a higher-frequency spatial mode, cf. Figs. 13c-e. The emergence of the higher-frequency mode is accompanied by a gradual decrease in both the population and wealth densities. At a slightly later time, the higherfrequency mode disappears and the spatial distribution again exhibits four distinct peaks but at different locations, so that the peaks and the troughs exchange places, cf. Figs. 13a and 13g. Further dynamics lead to the decay in both system components and eventually to system's extinction, see Figs. 13g-h. Figure 14 shows, for the same parameter values, the population and wealth densities vs time obtained at two fixed locations in space, i.e. at the boundary of the domain x = 0 and at x = (8/3)L, see Figs. 14a and c, respectively. It is readily seen that the change of the nearly-stationary dynamics (until approximately t = 17280) to much faster dynamics occurs when the local evolution of the system's variables takes them away from the vicinity of the saddle point (see Fig. 9).

Discussion and concluding remarks
In this paper, we have revisited factors and mechanisms affecting the spatial distribution of the human population. Such distributions often exhibit remarkable heterogeneity so that the population density in some areas (e.g. urban) can be much larger than in others (e.g. rural). In the modern perspective, there are many factors that contributes to this difference [27]. However, how the heterogeneity in the population distribution has developed over a longer, historical timescale is not always clear. A textbook explanation relates the emergence of densely populated areas to the heterogeneity of the environment. Centuries ago, the humankind was much more exposed to the forces of nature than it is nowadays. Areas with milder climatic conditions would have more likely been selected to establish a settlement. Recall also that the agriculture was the main driver of the economy, its efficiency being to a large extent dependant on the properties of the natural environment.
Convincing as it may sound, in this paper we endeavoured to challenge the above explanation. We first identified a few areas (selected from different parts of the world) where the environmental properties such as the average annual temperature, annual precipitation and the elevation do not show much variability in space over stretches of thousands of miles. We have revealed that, in spite of this apparent spatial homogeneity of the environment, the population density distribution over those areas is clearly heterogeneous -in fact, in all cases exhibiting a nearly-periodic spatial pattern.
We then ask the following question: can there exist another mechanism leading to the formation of a heterogeneous spatial distribution, even in an approximately homogeneous environment? We hypothesize that an appropriate mechanism can result from the nonlinear interplay between the human population and the resources that support its growth, in the manner of resource-consumer interaction [39]. Indeed, it is well known that, in a system of two or more interacting components, a locally stable steady state can become unstable with respect to a spatially heterogeneous perturbation with a certain wave length: the phenomenon called the Turing instability [51]. As a result, a spatially periodic pattern can arise [37].
Pattern formation due to the Turing instability is well known for chemical and biological systems [21,37]. For demographical systems, however, to the best of our knowledge the Turing instability has never been considered. In this work we have found this phenomenon within the framework of the demographic-economic model (15)(16) proposed in our earlier work [55,59]. Conditions of the Turing instability in terms of the model parameters are found analytically (see Section 5). The spatiotemporal dynamics of the system resulting from the instability, including the development of the periodic spatial pattern, was considered in computer simulations (Section 6).
Interestingly, our simulations reveal that the emerging periodical pattern is not the large-time asymptotics of the system, as it often happens in the case of Turing patterns, but the long-term transient dynamics [38]. After the pattern emerge (at time t ∼ 10 2 ), it remains almost unchanged over long time, on the order of t = 10 4 (for parameters of Fig. 13, until t ≈ 17200). Eventually, this quasi-stationary regime turns into a fast spatiotempotal dynamics where the emergence of a higher-frequency spatial mode is followed by the system collapse. Note that for the parameter values outside of the Turing instability range the spatial system can persist at its positive steady state indefinitely. Therefore, ultimately, the Turing instability drives the system to extinction.
Our study leaves a few open questions. Firstly, we mention that, although our theoretical findings (such as the formation of the Turing's spatial patterns) are in a qualitative agreement with the real-world data (see Section 2), a direct comparison between theory and data is hardly possible at this stage of research. Such comparison would required a sufficiently accurate estimate of the value of all model's parameters. This is a challenging and tedious task and will become a focus of a separate study. Secondly, our model is quite schematic. One question that arises here is what can be the effect of different ethnicity and/or race on the population dynamics. Different ethnicity can give rise to a different culture and that can affect the ways how the resources are consumed and the wealth is generated and distributed. Although we do not expect that it will change the system's properties completely -the existence of the nearly-periodic spatial pattern in countries as different as Canada, Australia and Mongolia points out at the universality and robustness of this phenomenon -yet this issue should be addressed more carefully, for instance by using more realistic models.
Another factor that was completely disregarded in our present study is the way how the given country's population has originally emerged. For instance, Canadas population mostly emerged as a result of migration from Europe, eventually propagating across Canada from East to West. In Australia, Europeans first settled in and around bays and then diffused over the continent. Finally, Mongolias population consists of mainly indigenous population. In mathematical terms, the different ways of how different countries came to be populated may correspond to different initial conditions; that, in its turn, may affect the selection of the emerging spatial mode.
Finally, the effect of the system's spatial dimensionality remains an open question. Recall that our spatially explicit model only include one spatial dimension. Arguably, it was appropriate for considering the population dynamics along a narrow stripe, as was the case identified in our realworld examples. Yet a mode general study on pattern formation in the demographic-economic system should consider it in the more realistic case of two spatial dimensions. That may reveal more complicated patterns and more complicated dynamics. That should be a focus of future research.