V Volpert
email: volpert@math.univ-lyon1.fr

S Petrovskii

A Zincenko

Interaction of human migration and wealth distribution

Keywords:

Dynamics of human populations depends on various economical and social factors. Their migration is partially determined by the economical conditions and it can also influence these conditions. This work is devoted to the analysis of the interaction of human migration and wealth distribution. The model consists of a system of equations for the population density and for the wealth distribution with conventional diffusion terms and with cross diffusion terms describing human migration determined by the wealth gradient and wealth flux determined by human migration. Wealth production and consumption depend on the population density while the natality and mortality rates depend on the level of wealth. In the absence of cross diffusion terms, dynamics of solutions is described by travelling wave solutions of the corresponding reaction-diffusion systems of equations. We show persistence of such solutions for sufficiently small cross diffusion coefficients. This result is based on the perturbation methods and on the spectral properties of the linearized operators.

Introduction

Dynamics of human population has long been a focus of research and controversy. Early insights into its uncontrolled growth resulted in a rather negative outlook and even disastrous conclusions because of the finiteness of available resources [START_REF] Malthus | An essay on the Principal of Population[END_REF] (see also Section 1.1 in [START_REF] Murray | Mathematical Biology I. An Introduction[END_REF]). Although nowadays this gloomy forecast is not widely shared, good understanding of the factors controlling the growth and proliferation of human population is required in order to avoid crises, poverty and social tensions.

Due to the effect of cultural, religious, political and other phenomena specific for human societies, e.g. see [START_REF] Turchin | War, space, and the evolution of Old World complex societies[END_REF], the population dynamics of humans is apparently much more complicated than the dynamics of animal species. Mathematical models of human dynamics allowing for the specifics and complexity of social and cultural interactions are at their infancy. There are, however, factors and processes that are important for all living beings, and one of them is movement. As well as animal species, the dynamics of human population has a distinct spatial aspect. In particular, human migration is a ubiquitous phenomenon that has been shaping countries and societies over centuries of human history [START_REF] Turney | Catastrophic early Holocene sea level rise, human migration and the Neolithic transition in Europe[END_REF][START_REF] Graf | Hunter-gatherer dispersals in the mammoth-steppe: technological provisioning and land-use in the Enisei River valley, south-central Siberia[END_REF] and that takes place at various scales ranging from individuals to nations [START_REF] Kupiszewski | Extra-union migration: the east-west perspective[END_REF][START_REF] Rogers | Introduction to Multiregional Mathematical Demography[END_REF]. The direction of migration (or, more generally, population flow) is often determined by the availability of resources and/or by the quality of living. 1 That can be thought about as the effect of the resource gradient. But the effect of resource gradient is known to be a factor that affects invasion (or, more generally, dispersal) of animal and plant species too, e.g. see [START_REF] Lortie | Effects of a directional abiotic gradient on plant community dynamics and invasion in a coastal dune system[END_REF][START_REF] Li | Success, failure, and spreading speeds for invasions on spatial gradients[END_REF]. Therefore, mathematical approaches that were developed to describe spatiotemporal dynamics of ecological populations [START_REF] Murray | Mathematical Biology I. An Introduction[END_REF][START_REF] Malchow | Spatiotemporal Patterns in Ecology and Epidemiology: Theory, Models, Simulations[END_REF] can be expected to be a reasonable starting point to describe the human dynamic too.

Proliferation of biological populations in space occurs due to an interplay between their reproduction and the movement of individuals [START_REF] Turchin | Quantitative Analysis of Movement[END_REF][START_REF] Hengeveld | Dynamics of Biological Invasions[END_REF][START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF]. The density p(x, t) of the population which depends on the space x and on time t can be calculated from a mathematical model. Nowadays, a broad variety of models can be used for this purpose [START_REF] Lewis | The Mathematics Behind Biological Invasions[END_REF] depending on the properties of the environment, e.g. uniform or fragmented [START_REF] Mistro | Spatiotemporal complexity of biological invasion in a space-and timediscrete predator-prey system with the strong Allee effect[END_REF], and on the type of individual movement, e.g. ordinary diffusion (Brownian motion) or anomalous diffusion [START_REF] Mendez | Stochastic Foundations in Movement Ecology[END_REF]. One commonly used type of model is given by the reaction-diffusion equations [START_REF] Volpert | Reaction-diffusion waves in biology[END_REF]. In the case where only a single species is considered explicitly, the model consists of a single partial differential equation:

∂p ∂t = D ∂ 2 p ∂x 2 + f (p), (1.1)
where the diffusion term corresponds to random movement of the individuals (ordinary diffusion), and the reaction term represents their natality and mortality. This equation with the function f (p) = p(1 -p) was introduced in [START_REF] Fisher | The wave of advance of advantageous genes[END_REF] in order to describe propagation of the dominant gene. In a more general ecological context, f (p) = αp(K -p) where α is the per capita population growth rate and K is the carrying capacity of the environment [START_REF] Murray | Mathematical Biology I. An Introduction[END_REF].

A remarkable property of Eq. (1.1) is that, if considered in an unbounded uniform space, for a broad class of initial conditions its large-time asymptotical solution is given by a travelling wave, that is by the solution p(x, t) = P (x -ct), where c = (2Dα) 1/2 is the speed of propagation [START_REF] Kolmogorov | A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem[END_REF]. Existence, stability and the speed of propagation of travelling waves are studied in detail for the logistic function f (p) and for more general nonlinearities (see [START_REF] Volpert | Reaction-diffusion waves in biology[END_REF][START_REF] Volpert | Traveling Wave Solutions of Parabolic Systems[END_REF] and the references therein).

Considering Eq. (1.1) in the context of population dynamics, it therefore predicts that the propagation of the population is described by a travelling wave with the constant speed c determined by the coefficients of the equation. The properties of the propagation becomes more complicated if other species (e.g. a predator or pathogen or competitor) or components are included (e.g. resource); the model is then given by a systems of PDEs. Well known examples are given by a system of competing species and by the prey-predator system where a simple travelling population front may turn into a sequence of fronts [START_REF] Petrovskii | Diffusive waves, dynamical stabilization and spatio-temporal chaos in a community of three competitive species[END_REF] or even a more complicated pattern of spread [START_REF] Volpert | Reaction-diffusion waves in biology[END_REF][START_REF] Morozov | Spatiotemporal complexity of the patchy invasion in a predator-prey system with the Allee effect[END_REF].

In the next section we will introduce the model of economical migration studied in this work. It consists of a system of two equations with cross diffusion terms. In the case of zero cross diffusion coefficients, it reduces to a conventional reaction-diffusion system of equations. We will study the case of small cross diffusion coefficients and will show that migration of the population can be described by travelling wave solutions.

Model of economical migration

In order to study economical migrations, we will complement the baseline model of biological invasions (1.1) with an equation for the resource or "wealth". Let p denote the population density and u the wealth distribution. Both of them are supposed to depend on space and time. We will consider the following system of equations:

∂u ∂t = D 1 ∂ 2 u ∂x 2 + D 2 ∂ ∂x  u ∂p ∂x  + W (u, p) -S(u, p), ( 2.2
)

∂p ∂t = D 3 ∂ 2 p ∂x 2 -D 4 ∂ ∂x  p ∂u ∂x  + αp n (K(u) -p) -σ(u)p. (2.3)
Here the diffusion term in the first equation describes local wealth redistribution due to the economical activity or taxes, the second diffusion term describes wealth transport due the displacement of individual, W and S are the wealth production and consumption. These functions will be specified below. The first diffusion term in the second equation describes random motion of the population, and the second diffusion term shows the migration of the population in the direction of wealth gradient. The second term in the right-hand side of this equation describes reproduction of the population where the power n is introduced in order to attain generality and to account for collective effects. The special case of n = 1 corresponds to the standard logistic population growth [START_REF] Murray | Mathematical Biology I. An Introduction[END_REF]. The last term in Eq. (2.3) describes mortality. In a general case, carrying capacity K and mortality rate σ depend on wealth. Let us note that cross diffusion terms resemble the Soret and Dufour effects in heat and mass transport, the second diffusion term in Eq. ( 2.3) is similar to that in chemotaxis. Production of wealth is often considered in the form of Cobb-Douglas production function:

W (L, Q, M ) = bL α Q β M γ , ( 2.4)
where L i represents the labour input, Q i the capital, machinery, transportation, infrastructure, M i natural resources and materials [START_REF]The digital economist[END_REF]; a positive coefficient b represents a measure of technology, α, β and γ are positive constants. The condition that doubling the inputs leads to doubling the output signifies that α + β + γ = 1. This model satisfies the law of diminishing marginal productivity [START_REF] Groth | Lecture notes in macroeconomics[END_REF]. In some models, this sum is considered different from 1. Assuming that M is constant, capital Q is proportional to wealth and labour L is proportional to the population density, we get the wealth production rate in the form W (u, p) = f (u)g(p). We will consider the functions f (u) and g(p) in a more general form. Typically these are growing functions with saturation.

Next, we specify the rate of wealth consumption. At the level of individual consumption c, we will consider Keynes consumption function, c = r +sy, where y is income, r and s are some positive coefficients. Assuming that average income is proportional to the wealth, we will consider the rate of wealth consumption in the form S(u, p) = au + (r + su)p, where the first term corresponds to its consumption independent of the population level and the second term depending on it.

Let us now estimate the dependence of carrying capacity K(u) on wealth u. Suppose that the number of children N in the family is proportional to the family income R and inversely proportional to the child related expenses E, N ∼ R/E. Then there are two limiting cases. In the first one, E is fixed, and R is proportional to the wealth u. In this case, K(u) is a growing function of u. In the second case, R is still proportional to u, but E also depends on u and it grows faster than R. In particular, E dependence on u can be related to education and various commodities required for children in modern society. In this second case, K(u) is a decreasing function of u. Both patterns are known. When modern technologies come to undeveloped countries, they lead to a rapid population growth. On the other hand, it is also known that natality rate decreases with growth of the education level. Finally, the mortality rate decreases as a function of wealth, σ(u) = σ 0 -σ 1 u/ (1 + u).

In this work we study the existence of travelling waves for system (2.2), (2.3) in the case of small cross diffusion coefficients. If D 2 = D 4 = 0, then the wave existence follows from the properties of the monotone systems [START_REF] Volpert | Elliptic Partial Differential Equations[END_REF][START_REF] Volpert | Traveling Wave Solutions of Parabolic Systems[END_REF]. In order to study the case where these parameters are positive and sufficiently small, we apply the perturbation technique based on the implicit function theorem. We begin with the spectral properties of the corresponding linear operators and establish the conditions under which they satisfy the Fredholm property (Section 3). It allows us to formulate the solvability conditions and prove the invertibility of linear operators on properly chosen subspaces. The invertibility is used in the implicit function theorem. In Section 4 we employ it to prove the wave existence with a variable speed in the bistable case (Theorem 4.1) and with a fixed speed in the monostable case (Theorem 4.2). These results are applied to study the migration model.

Spectral properties of linear operators

We consider the system (2.2), (2.3)

∂u ∂t = D 1 ∂ 2 u ∂x 2 + D 2 ∂ ∂x  u ∂p ∂x  + F (u, p), ( 3.1
)

∂p ∂t = D 3 ∂ 2 p ∂x 2 -D 4 ∂ ∂x  p ∂u ∂x  + G(u, p) (3.2)
with more general functions F (u, p) and G(u, p) in view of other possible applications. We will consider travelling wave solutions of this systems, that is solutions of the form u(x, t) = U (x -ct), p(x, t) = P (x -ct), which satisfy the system of equations

D 1 U ′′ + D 2 (U P ′ ) ′ + cU ′ + F (U, P ) = 0, (3.3
)

D 3 P ′′ -D 4 (P U ′ ) ′ + cP ′ + G(U, P ) = 0 (3.4)
considered on the whole axis. Here prime denotes the derivative with respect to the variable z = x -ct.

Consider the operator A corresponding to the left-hand side of system (3.3), (3.4):

A  U P  =  D 1 U ′′ + D 2 (U P ′ ) ′ + cU ′ + F (U, P ) D 3 P ′′ -D 4 (P U ′ ) ′ + cP ′ + G(U, P )
and the operator L linearized about some functions U 0 , P 0 :

L  U P  =  a 11 U ′′ + a 12 P ′′ + b 11 U ′ + b 12 P ′ + c 11 U + c 12 P a 21 U ′′ + a 22 P ′′ + b 21 U ′ + b 22 P ′ + c 21 U + c 22 P ,
where

a 11 = D 1 , a 12 = D 2 U 0 (z), a 21 = -D 4 P 0 (z), a 22 = D 3 , (3.5) b 11 = c + D 2 P ′ 0 (z), b 12 = D 2 U ′ 0 (z), b 21 = -D 4 P ′ 0 (z), b 22 = c -D 4 U ′ 0 (z), (3.6
)

c 11 = F ′ U (U 0 , P 0 ) + D 2 P ′′ 0 , c 12 = F ′ P (U 0 , P 0 ), c 21 = G ′ U (U 0 , P 0 ), c 22 = G ′ P (U 0 , P 0 ) -D 4 U ′′ 0 . (3.7)
We suppose that the diffusion coefficients D i , i = 1, 2, 3, 4 are positive, the functions F (U, P ) and G(U, P ) have continuous first partial derivatives with respect to U and P , and these derivatives satisfy the Lipschitz condition in every bounded set in R 2 . Moreover, we assume that U 0 , P 0 ∈ C 2+α (R) and U 0 (z), P 0 (z) > 0 for all z ∈ R. Then the coefficients of the operator L belong to the space C α (R). We consider the operator L as acting from the space

E 1 = C 2+α (R) into the space E 2 = C α (R).
Let us recall that the operator L satisfies the Fredholm property if its kernel has a finite dimension, its image is close and the codimension of the image has also a finite dimension. The latter means that the nonhomogeneous equation Lw = f is solvable if and only if a finite number of solvability conditions φ i (f ) = 0, i = 1, . . . , N are satisfied. Here w = (U, P ), φ i are some functionals from the dual space E * 2 . In order to study the Fredholm property of this operator, we suppose that the functions U 0 and P 0 have limits at infinity and introduce the limiting operators

L ±  U P  =  a ± 11 U ′′ + a ± 12 P ′′ + b ± 11 U ′ + b ± 12 P ′ + c ± 11 U + c ± 12 P a ± 21 U ′′ + a ± 22 P ′′ + b ± 21 U ′ + b ± 22 P ′ + c ± 21 U + c ± 22 P,
where the superscripts in the coefficients signify their limiting values at ±∞.

Proposition 3.1. Suppose that the matrix A(z) = (a ij (z)
) is such that det A(z) ≥ a 0 for all z ∈ R and for some positive number a 0 . Then the operator L is normally solvable with a finite dimensional kernel if and only if the limiting equations

L ± w = 0 (3.8)
have only zero solutions in E 1 .

The proof of the proposition follows from more general results on elliptic operators in unbounded domains [25, Chapter 4, Theorems 3.1, 3.2]. Let us note that by virtue of the condition on the determinant of the matrix A, the operator L is elliptic in the Douglis-Nirenberg sense. Taking into account the explicit form (3.5) of the coefficients, we get

det A = D 1 D 3 + D 2 D 4 U 0 P 0 ≥ D 1 D 3 > 0.
Therefore the condition of the proposition is satisfied.

Since the limiting operators have constant coefficients, then we can apply the Fourier transform to equations (3.8). These equations have only zero solutions if and only if

det  -ξ 2 A ± + iξB ± + C ±  ̸ = 0 ∀ξ ∈ R, (3.9)
where

A ± =  a ± 11 a ± 12 a ± 21 a ± 22  , B ± =  b ± 11 b ± 12 b ± 21 b ± 22  , C ± =  c ± 11 c ± 12 c ± 21 c ± 22  . Set M ± (ξ) = -ξ 2 A ± + iξB ± + C ± .
Taking into account the explicit form of the coefficients (3.5)-(3.7), we get

M ± (ξ) =  -ξ 2 D 1 + s ± 11 , -ξ 2 D 2 U ± 0 + s ± 12 ξ 2 D 4 P ± 0 + s ± 21 , -ξ 2 D 3 + s ± 22  + ciξI,
where I is the identity matrix and

s ± 11 = F ′ U (U ± 0 , P ± 0 ), s ± 12 = F ′ P (U ± 0 , P ± 0 ), s ± 21 = G ′ U (U ± 0 , P ± 0 ), s ± 22 = G ′ P (U ± 0 , P ± 0 ).
In what follows, the superscripts ± will be omitted for brevity. Consider first the case c = 0. Then det M (ξ) = αξ 4 + βξ 2 + γ,

where α = D 1 D 3 + D 2 D 4 U 0 P 0 , β = -D 1 s 22 -D 3 s 11 + D 2 U 0 s 21 -
or β 2 ≥ 4αγ and -β +  β 2 -4αγ < 0 (3.11)
is satisfied and for c ̸ = 0 conditions γ ̸ = 0 and det M (ξ 0 ) ̸ = 0, (3.12)

where ξ 2 0 = (s 11 + s 22 )/(D 1 + D 3 ), are satisfied.
The proof of this theorem is given in Appendix A.

Remark 3.3. Let us recall that the essential spectrum of the operator L is the set of all complex numbers λ for which the operator L -λ does not satisfy the Fredholm property. According to the results presented above, it is determined by the equality

det(M (ξ) -λI) = 0, ξ ∈ R. (3.13)
If the essential spectrum of the operator L lies completely in the left-half plane, then the index of the operator equals 0. Indeed, the operator L -λ is invertible for large positive λ [START_REF] Volpert | Elliptic Partial Differential Equations[END_REF], and the index is constant in each connected set at the complex plane where the operator satisfies the Fredholm property.

If det M (ξ) = 0 for some real ξ, then the essential spectrum of the operator L crosses the origin, and the operator does not satisfy the Fredholm property. This is also related to the instability of solutions by means of spatial structures emerging at infinity. Suppose that the matrix S = (s ij ) has both eigenvalues with the negative real part. Then γ > 0 and s 11 + s 22 < 0. This condition is not sufficient to provide that det M (ξ) ̸ = 0. Indeed, if the coefficient β is negative and sufficiently large in the absolute value, then det M (ξ) can vanish. In the case of diagonal diffusion matrix (D 2 = D 4 = 0), if one of the coefficients s 11 and s 22 is positive and the corresponding diffusion coefficient is sufficiently large, then det M (ξ) = 0 for some ξ. This is related to the Turing instability. If both coefficients s 11 and s 22 are negative, then the equality det M (ξ) = 0 can be achieved due to the cross diffusion coefficients if s 12 > 0 or s 21 < 0.

We also observe that conditions (3.10), (3.11) are not satisfied for γ ≤ 0. Indeed, this is obvious for condition (3.10). As to condition (3.10), since α > 0, then -β +  β 2 -4αγ ≥ 0.

Small cross diffusion coefficients

Perturbations of waves

In the case of small cross diffusion coefficients we set D 2 = ϵd 2 , D 4 = ϵd 4 , where ϵ, d 2 and d 4 are positive numbers. System (3.3), (3.4) writes

D 1 U ′′ + ϵd 2 (U P ′ ) ′ + cU ′ + F (U, P ) = 0, ( 4.1)
D 3 P ′′ -ϵd 4 (P U ′ ) ′ + cP ′ + G(U, P ) = 0. (4.2)
We will prove the existence of solutions of this system for sufficiently small ϵ assuming that the solution exists for ϵ = 0. We will consider below applications to the migration model where the existence for ϵ = 0 can be verified.

Suppose that F (U ± , P ± ) = G(U ± , P ± ) = 0 for some (U ± , P ± ) and set

S ± =  F ′ U (U ± , P ± ) F ′ P (U ± , P ± ) G ′ U (U ± , P ± ) G ′ P (U ± , P ± )  .
L 0  U P  =  D 1 U ′′ + c 0 U ′ + F ′ U (U ± , P ± )U + F ′ P (U ± , P ± )P D 3 P ′′ + c 0 P ′ + G ′ U (U ± , P ± )U + G ′ P (U ± , P ± )P
is simple, then for all ϵ sufficiently small system (4.1), (4.2) has a solution for some c = c ϵ .

Proof. We consider the operator

A ϵ  U P  =  D 1 U ′′ + ϵd 2 (U P ′ ) ′ + cU ′ + F (U, P ) = 0, D 3 P ′′ -ϵd 4 (P U ′ ) ′ + cP ′ + G(U, P ) = 0
corresponding to system (4.1), (4.2). It acts from the space E 1 × R into the space E 2 , where

E 1 = C 2+α (R), E 2 = C α (R) (U, P ∈ E 1 , c ∈ R).
According to the conditions of the theorem, equation A 0 (w) = 0 has a solution w = w 0 . Here w = (U, P ), w 0 = (U 0 , P 0 ). We will use the implicit function theorem in order to prove the existence of solutions of the equation A ϵ (w) = 0 for all sufficiently small ϵ > 0. It can be easily verified that the operator A ′ ϵ (w 0 ) is continuous with respect to (ϵ, w 0 ) (see, e.g., [START_REF] Apreutesei | Travelling waves for integro-differential equations in population dynamics[END_REF]). It remains to verify the invertibility of the operator L 0 .

Let us begin with the spectral properties of the operator L 0 . The essential spectrum of this operator is given by the equality (3.13) where

det(M ± (ξ) -λI) = α ± ξ 4 + (β ± -c 2 0 )ξ 2 + γ ± + λ 2 + (c 0 iξ -λ)  s ± 11 + s ± 22 -(D 1 + D 3 )ξ 2  .
Since the matrices S ± have eigenvalues with negative real parts, then γ ± = det S ± > 0, s ± 11 + s ± 22 < 0. It can be directly verified that the essential spectrum λ(ξ) of the operator L 0 determined as a solution of Eq. (3.13) lies in the left-half plane of the complex plane for all real ξ. Therefore the operator L 0 satisfies the Fredholm property and it has the zero index (see Remark 3.3). Differentiating Eqs. (4.1), (4.2), we obtain that the operator L 0 has a zero eigenvalue with the eigenfunction (U ′ 0 , P ′ 0 ). Therefore it is not invertible and we cannot apply the implicit function theorem directly. We will show that it is applicable on a properly chosen subspace.

Equation L 0 w = f is solvable if and only if  ∞ -∞ f (z)w * (z)dz = 0,
where w * is the eigenfunction corresponding to the zero eigenvalue of the formally adjoint operator L * .

Let us consider the subspace E 0 1 of the space

E 1 = C 2+α (R) which consists of functions from E 1 for which  ∞ -∞ f (z)w ′ 0 (z)dz = 0,
where w ′ 0 = (U ′ 0 , P ′ 0 ). Since the eigenfunction w ′ 0 is exponentially decaying at infinity, then this integral is well defined. We consider the operator A ϵ as acting from

E 0 1 × R (U, P ∈ E 0 1 , c ∈ R) into E 2 = C α (R)
. The constant c is unknown and it should be taken into account in the linearized operator:

A ′ ϵ (w 0 )w = L 0 w + cw ′ 0 ,
where w 0 = (U 0 , P 0 ), w = (U, P ).

From the condition that the zero eigenvalue of the operator L 0 is simple, it follows that

 ∞ -∞ w ′ 0 (z)w * (z)dz ̸ = 0. (4.3)
Indeed, if this integral equals zero, then the equation L 0 w = w ′ 0 has a solution w 1 . Therefore L 0 w 1 ̸ = 0, L 2 w 1 = 0, and the zero eigenvalue is not simple.

Hence equation A ′

ϵ (w 0 )w = f is solvable for any f ∈ E 2 . Indeed, if we write it in the form

L 0 w = f -cw ′ 0 ,
then we determine c from the equality

 ∞ -∞ f (z)w * (z)dz = c  ∞ -∞ w ′ 0 (z)w * (z)dz (cf. (4.
3)), and the solvability condition is satisfied.

Thus, the operator

A ′ ϵ (w 0 ) : E 0 1 × R → E 2 is invertible.
The assertion of the theorem follows from the implicit function theorem. Theorem 4.2. Suppose that system (4.1), (4.2) has a solution (U 0 (z), P 0 (z)) with the limits U 0 (±∞) = U ± , P 0 (±∞) = P ± for ϵ = 0 and some c = c 0 ̸ = 0. Suppose, next, the matrix S + has one positive and one negative eigenvalue, and the matrix S -has both eigenvalues in the left-half plane. If the dimension of the kernel of the operator L 0 equals 1 and the codimension of its image equals 0, then for all ϵ sufficiently small system (4.1), (4.2) has a solution for the same value c = c 0 .

Proof. Similar to the proof of the previous theorem we consider the operator

A ϵ : E 1 × R → E 2 .
According to the assumption of the theorem, equation A 0 (w) = 0 has a solution w 0 (w = (U, P ), w 0 = (U 0 , P 0 )). We will use the implicit function theorem in order to prove that equation A ϵ (w) = 0 has a solution for all ϵ > 0 sufficiently small.

From the conditions of the theorem it follows that the linearized operator L 0 = A ′ 0 (w 0 ) satisfies the Fredholm property, the dimension of its kernel equals 1, and its image coincides with the whole space E 2 . Let E 0 1 be the subspace of the space E 1 = C 2+α (R) which consists of functions from E 1 orthogonal in L 2 to the eigenfunction (U ′ 0 , P ′ 0 ). We consider the operator A ϵ as acting from E 0 1 into E 2 = C α (R). Since the operator L 0 : E 0 1 → E 2 is invertible, then we can apply the implicit function theorem. Therefore equation A ϵ (w) = 0 has a solution w ∈ E 0 1 for all ϵ > 0 sufficiently small. Remark 4.3. Existence of waves is proved in the bistable case (Theorem 4.1) and in the monostable case (Theorem 4.2). In the former, the wave speed depends on ϵ, while in the latter, the wave speed is fixed. This difference is determined by the index of the linearized operators. It is 0 in the bistable case and 1 in the monostable case. We will consider the applications of these theorems in the next section.

Applications to the migration model

In order to apply the results of the previous section we begin with the case of zero cross diffusion coefficients:

D 1 U ′′ + cU ′ + F (U, P ) = 0, (4.4
)

D 3 P ′′ + cP ′ + G(U, P ) = 0, (4.5)
where

F (U, P ) = f (U )g(P ) -au, G(U, P ) = αP n (K(u) -P ) -σP.
We assume here for simplicity that a and σ are constant. More general case presented in the introduction can also be considered. We have also introduced the power n of the population density in the reproduction term. If n = 0 (and eventually σ = 0), then we neglect the reproduction and mortality and suppose that the evolution of the population is determined by the carrying capacity K(u). In the case n = 1 we obtain conventional logistic equation where the reproduction rate depends only on female density, if n = 2 then it depends on the product of male and female densities assumed to be equal to each other. In what follows we suppose that f ′ (U ) > 0, g ′ (P ) > 0, K ′ (U ) > 0, that is wealth production increases with wealth and population density, carrying capacity increases as a function of wealth. Therefore

∂F ∂P > 0, ∂G ∂U > 0, (4.6)
system (4.4), (4.5) is monotone, and we can apply the results on wave existence [START_REF] Volpert | Traveling Wave Solutions of Parabolic Systems[END_REF].

Bistable case without reproduction. For simplicity, we begin with the case n = 0 and set α = 1, σ = 0. Note that this case is apparently not entirely biologically realistic, because it assumes that the local population growth rate can be positive for P = 0 as G(U, 0) ̸ = 0 (which can as well mean that the system is not closed and there are other sources of population migration that we do not discuss here). However, this case is a convenient example to demonstrate the mathematical technique and to reveal the principle properties of the system.

The zeros of the nonlinearities can be found from the equations f (U )g(P ) -aU = 0, K(U ) -P = 0.

We suppose that f (U ) is an S-shape function qualitatively shown in Fig. 1. The functions g(P ) and K(U ) have a similar behaviour with possibly different inflection points. Set h(U ) = f (U )g(K(U )). Equation h(U ) = aU can have different numbers of solutions. Two possible situations are shown in Fig. 1. In the first case, there are three solutions, the points (U ± , P ± ) are stable in the sense that the corresponding matrices (Theorem 4.1) have eigenvalues in the left-half plane. The intermediate point is unstable, the corresponding matrix has a positive eigenvalue. There exists a unique solution (U 0 (z), P 0 (z)) of system (4.4), (4.5). All conditions of Theorem 4.1 are satisfied, so there exists a solution for small cross diffusion coefficients.

In the second case, there are five solutions. Besides the points (U ± , P ± ), there is one more stable point (U * , P * = K(U * )). In this case, system (4.4), (4.5) has a solution (U 1 (z), P 1 (z)) with the limits

U 1 (-∞) = U -, P 1 (-∞) = P -, U 1 (+∞) = U * , P 1 (+∞) = P * (4.7)
and another solution (U 2 (z), P 2 (z)) with the limits

U 1 (-∞) = U * , P 1 (-∞) = P * , U 1 (+∞) = U + , P 1 (+∞) = P + . (4.8)
The corresponding values of the wave speed are c 1 and c 2 , respectively. If c 1 > c 2 , then there exists a solution (U 0 (z), P 0 (z)) with the limits (U ± , P ± ). This solution describes the dynamics of the original system (2.2), (2.3) (without the cross-diffusion terms) in the large time limit. Similar solutions exist in the case of small nonzero cross diffusion coefficients. In case c 1 ≤ c 2 , such solution does not exist [START_REF] Volpert | Wave trains described by monotone parabolic systems[END_REF][START_REF] Volpert | Elliptic Partial Differential Equations[END_REF]. For intermediate time, the dynamics of the system can more generally be given by a system of two travelling front; see Fig. 2.

Bistable case with reproduction. In the case n = 2 we will suppose for simplicity that g(P ) and K(u) are linear functions, g(P ) = P, K(U ) = U , and α = 1. This assumption allows us to get the analytical expressions for the zero line of the functions F and G: P = aU f (U ) , U = P + σ P (P ̸ = 0) and P = 0 (Fig. 3). The number of stationary points of the kinetic system of equations depends on the parameters. Two examples are shown in Fig. 3. There are two stable stationary points in the left figure, (U + , P + ) = (0, 0) and (U -, P -) > 0. Other points are unstable. Therefore there exists a unique (up to a shift) monotonically decreasing solution of system (4.4), (4.5) with the limits (U ± , P ± ) at infinity. There are three stable stationary points in the example in the right figure. Denote the intermediate point by (U * , P * ). Then there are waves with the limits (4.7), (4.8). Denote their speeds by c 1 and c 2 , respectively. As above in the case n = 0, the wave with the limits (U ± , P ± ) exists if and only if c 1 > c 2 . In all these cases, the waves also exist for sufficiently small cross diffusion coefficients.

dU dt = F (U, P ), dP dt = G(U, P )
Monostable case with reproduction. In the case n = 1, the zero lines of the functions F and G are determined by the system of equations:

f (U )g(P ) -aU = 0, P = K(U ) - σ α and P = 0.
The number and stability of the corresponding stationary points depend on the parameters and on the form of the functions f (U ), g(P ) and K(U ). We will not analyse here all possible cases. Let us note that the point with P = 0 is unstable if αK ′ (U 0 ) > σ, where U 0 is a solution of the equation F (U )g(0) = aU . Suppose that the point (U + , P + ) is unstable and the point (U -, P -) is stable in the sense that the corresponding matrix has a positive eigenvalue in the first case, and both of them are negative in the second case. Moreover, we assume that U -> U + , P -> P + , and there are no other stationary points in the rectangle U + ≤ U ≤ U -, P + ≤ P ≤ P -. Then monotonically decreasing solutions of system (4.4), (4.5) with the limits (U ± , P ± ) exist for all values of speed c ≥ c 0 , where c 0 is some positive number.

The index of the linearized operator L 0 equals 1 [START_REF] Collet | Computation of the index of linear elliptic operators in unbounded cylinders[END_REF]. Therefore, if the dimension of the kernel of the operator is 1, then the codimension of the image is 0, and we can apply Theorem 4.2. In order to determine the dimension of the kernel of the operator L 0 , let us recall that the derivative (U ′ 0 (z), P ′ 0 (z)) of the solution is the eigenfunction corresponding to the zero eigenvalue. Assuming that all solutions λ i ± of the equations det T ± (λ) = 0, where T ± (λ) = Dλ 2 + cλ + S ± , D is the diagonal matrix with the diagonal elements D 1 and D 3 , are different from each other, we get

 U ′ 0 (z) P ′ 0 (z)  ∼ 4  i=1 k i ± q i ± e λ i ± z , z → ±∞. (4.9)
Here q i ± are some vectors, k i ± are some constants. It can be shown that there are three negative values λ i + which determine behaviour of solutions at +∞ (Appendix B). Assuming for simplicity that they are different from each other, λ 3 + < λ 2 + < λ 1 + < 0, we can verify that the waves with the speed c > c 0 decay with the maximal exponent λ 1 + . In this case, the kernel of the linearized operator has dimension 1, Theorem 4.2 is applicable, and we prove the wave existence for small cross diffusion coefficients.

If c = c 0 , then solution (U 0 (z), P 0 (z)) of system (4.4), (4.5) decays with the exponent λ 2 + < λ 1 + . In this case we can introduce an exponential weight at ∞ with some weight µ, λ 2 + < µ < λ 1 + . The essential spectrum of the linearized operator L 0 considered in the weighted space is completely in the left-half plane, and its index equals zero. The derivative (U ′ 0 (z), P ′ 0 (z)) of the solution belongs to the weighted space. Hence the dimension of the kernel of the operator and the codimension of its image equal 1. Similar to the bistable case, the solution persists under small perturbations of the system but the speed c 0 can change.

Discussion

In this work we study the interaction of human population with wealth distribution. This interaction concerns the fluxes and the production terms for both variables. Conventional diffusion term in the equation for the population density describes random motion of the individuals while the second diffusion terms corresponds to the migration of human population along the gradient of wealth distribution. Similar terms in the equation for wealth distribution describe diffusive flux and wealth transport due to human migration. Wealth transfer from rich to poor regions occurs by redistribution mechanisms and by various types of economic activity (trade, investment, delocalization, and so on). Production terms depend on both variables. Natality and mortality rates in the equation for the population density depend on wealth, production and consumption of wealth depend on the population density.

In the case of zero cross diffusion coefficients we obtain a reaction-diffusion system of equations. It satisfies the conditions of monotonicity (positive off-diagonal elements in the Jacobian of production rates). Such systems satisfy the maximum and comparison principle allowing the investigation of the existence and stability of travelling waves.

In the bistable case, if there are only two stable stationary points, then there exists a unique travelling wave connecting them. It is globally asymptotically stable. If there are three stable stationary points, then there are two consecutive waves connecting them. The relation between their speeds determines the existence of the merged wave. We prove the wave existence in the case of small cross diffusion coefficients using the implicit function theorem. Its application is based on spectral properties of the linearized operators. The perturbation of the problem influences the wave speed.

In the monostable case with zero cross diffusion coefficients the waves exist for all values of the speed greater than or equal to the minimal speed. The index of the corresponding linear operator equals 1. We show that its kernel has the same dimension. Therefore the image of the operators coincides with the whole space. This property of the operator allows us to prove wave existence for a fixed wave speed. Its value does not change for small nonzero cross diffusion coefficients; see Theorem 4.2 and Corollary 7 (Appendix B).

Note that, in this paper, it was not our intention to provide an exhausting study of human migrations. Our goal is more modest, i.e. to introduce a new spatially explicit model of coupled population-wealth dynamics and to reveal some of its basic mathematical properties, thus laying the foundation for future research in this direction. Correspondingly, we have only studied in this work dynamics of solutions under some particular conditions on diffusion and production terms. We can expect new effects and more complicated dynamics for other values of parameters and production functions. In particular, cross diffusion terms can lead to the emergence of patterns (cf. Remark 3.3 in Section 3). These questions will be studied in the subsequent works.

Hence for k sufficiently large we have β 2 τ < 4α τ γ τ . Hence the operator L τ is normally solvable with a finite dimensional kernel for all τ ∈ [0, 1]. Its index is independent of τ .

A similar construction can be done for c = 0 with condition (A.2) and for c ̸ = 0. Let us note that during the deformation condition (A.2) can change to (A.1). Then we proceed as above. In order to consider behaviour at +∞ and -∞ independently of each other we can introduce variable coefficients D 1 and D 3 with different limits at infinities. -and λ 2 -. If c > c 0 , then they are different, and a positive vector q i -in (4.9) at -∞ corresponds to the minimal of them.

Proof. If D 1 = D 3 , then the solutions of the equation det T -(λ) = 0 can be found explicitly from the equations

D 1 λ 2 + cλ + µ 1,2 -= 0, where µ 1,2
-< 0 are the eigenvalues of the matrix S -. In this case there are two positive and two negative values of λ. When we change the value of D 3 , these solutions cannot cross the imaginary axis since det T ± (iξ) ̸ = 0 for any real ξ provided that the essential spectrum does not cross the origin. Hence there are two values of λ with positive real parts. By virtue of the existence of monotonically decreasing solution (U 0 (z), P 0 (z)) these solutions are real. They do not coincide because such solutions persist for a small change of c due to the assumption of the lemma.

It remains to verify that positive vector q i -in (4.9) corresponds to the minimal of these solutions. Since T -(λ)q -= 0, then q -> 0 is the eigenvector corresponding to the zero eigenvalue of the matrix T -(λ). Hence 0 is the maximal eigenvalue of this matrix, and another one is negative. The matrix T -(0) has two negative eigenvalues. Therefore, T -(λ) has one negative and one zero eigenvalue for the minimal solution of the equation det T -(λ) = 0, and one positive and one zero eigenvalue for the maximal solution.

Lemma 4. Equation det T + (λ) = 0 has three solutions with negative real parts. If they are real and different, then positive vector q i + in (4.9) at ∞ corresponds to one of the two maximal solutions. The proof of the lemma is similar to the previous one.

Proposition 5. Suppose that all solutions of the equations det T ± (λ) = 0 are real and different. If the positive vector q i + in (4.9) at ∞ corresponds to the maximal negative solution of the equation det T + (λ) = 0, then the kernel of the operator L 0 has dimension 1, and c > c 0 .

Proof. The vector-function w(z) = -(U ′ 0 (z), P ′ 0 (z)) is a positive eigenfunction corresponding to the zero eigenvalue of the operator L 0 . According to Lemmas 3 and 4, we have w(z) ∼ k 1 + q 1 + e λ 1 + + k 2 + q 2 + e λ 2 + + k 3 + q 3 + e λ 3 + , z → ∞,

w(z) ∼ k 1 -q 1 -e λ 1 -+ k 2 -q 2 -e λ 2 -, z → -∞,
where k 1 ± ̸ = 0, q 1 ± > 0, 0 > λ 1 + > λ 2 + > λ 3 + , 0 < λ 1 -< λ 2 -. Suppose that there is another bounded solution w * (z) of the equation L 0 v = 0. Then w * (z) ∼ c 1 + q 1 + e λ 1 + + c 2 + q 2 + e λ 2 + + c 3 + q 3 + e λ 3 + , z → ∞,

w * (z) ∼ c 1 -q 1 -e λ 1 -+ c 2 + q 2 -e λ 2 -, z → -∞
with some constants c i ± . Without loss of generality we can assume that c 1 + = 0. Otherwise we can take a linear combination of the functions w * (z) and w(z) which satisfies this condition. Therefore w * (x) decays faster at ∞ than w(z). We can also assume that the vector-function w * (z) has some positive values for at least one of its components and for some z. Indeed, otherwise we can multiply it by -1.

Consider the function τ w(z) with some positive number τ . Since the matrix S -has negative eigenvalues, then there exists such z = z 0 that from the inequality τ w(z 0 ) ≥ w * (z 0 ) it follows that τ w(z) ≥ w * (z) for all z ≤ z 0 [START_REF] Volpert | Traveling Wave Solutions of Parabolic Systems[END_REF]. The inequalities between the vectors are understood component-wise. On the other hand, τ w(z) > w * (z) for any τ > 0 and all z sufficiently large (possibly depending on τ ). Hence τ w(z) > w * (z), z ∈ R (B.1)

for τ large enough.

Denote by τ 0 infimum of all τ for which τ w(z) > w * (z) for z ≥ z 0 . Then τ 0 > 0 since w * (z) has positive values. Inequality (B.1) holds for such values of τ . Clearly, τ 0 w(z) ≥ w * (z) for z ≥ z 0 and τ 0 w(z 1 ) = w * (z 1 ) for some z 1 ≥ z 0 and for one of the components of the vector-function. Therefore, τ 0 w(z) ≥ w * (z) for all z ∈ R. Thus, we obtain a contradiction with the positiveness theorem (or maximum principle) at z = z 1 .

We proved that w(z) is a unique up to a factor solution of equation L 0 v = 0. Since the index of this operator equals 1, then the image of the operator coincides with the whole space E 2 . The operator L 0 is invertible as acting from the subspace E 0 1 of functions orthogonal to w 0 into E 2 . From the implicit function theorem it follows that solution of system (4.4), (4.5) persists under a small variation of c. Hence c > c 0 . Corollary 6. A monotonically decreasing solution of system (4.4), (4.5), which decays at ∞ with the exponent λ 1 + < 0, (the maximal negative solution of the equation det T + (λ) = 0) is unique up to a shift and it corresponds to a speed c > c 0 . Solution with the minimal speed c 0 decays with the exponent λ 2 + < λ 1 + , which is the second solution of this equation.

Theorem 4 . 1 .

 41 Suppose that system (4.1), (4.2) has a solution (U 0 (z), P 0 (z)) with the limits U 0 (±∞) = U ± , P 0 (±∞) = P ± for ϵ = 0 and some c = c 0 , and one of the conditions (3.10), (3.11) is satisfied if c 0 = 0. If the matrices S + and S -have all eigenvalues in the left-half plane and the zero eigenvalue of the operator

Fig. 1 .

 1 Fig. 1. Qualitative form of the functions f (U ) and h(U ). Solutions of the equation h(U ) = U determine the stationary points and their stability.

Fig. 2 .

 2 Fig. 2. A sketch of the two wave system arising in the bistable case if c 1 ≤ c 2 (solid line corresponds to the P component of the solution, dashed line to the U component).

Fig. 3 .

 3 Fig. 3. Two specific examples of zero lines of the functions F and G. There are two stable points in the left figure (marked by "s") and three of them in the right figure.

Appendix B Lemma 3 .

 3 Equation det T -(λ) = 0 has two positive solutions λ 1

Corollary 7 .

 7 System (4.1), (4.2) has solutions for all c > c 0 and ϵ > 0 sufficiently small.

  Theorem 3.2. The operator L with the coefficients (3.5)-(3.7) satisfies the Fredholm property if and only if for c = 0 one of the conditions

	β 2 < 4αγ,	(3.10)

D 4 P 0 s 12 , γ = s 11 s 22 -s 12 s 21 .

Considering, for instance, the recent surge in the migration to Europe, whilst some of the migrants were apparently fleeing war or genocide, many migrants openly stated that they were in search of economic benefits and more comfortable life.

Acknowledgements

The first author was partially supported by the Ministry of Education and Science of the Russian Federation (the agreement number 02.a03.21.0008) and by the program PICS CNRS 6583 Matbio.

Appendix A

We will present here the proof of Theorem 3.2. We begin with auxiliary lemmas. Lemma 1. If c = 0, then the determinant of the matrix M (ξ) is different from 0 for all real ξ if and only if one of the following two conditions is satisfied:

or

Indeed, we consider det M (ξ) as a second-order polynomial with respect to ξ 2 . If condition (A.1) is satisfied, then it does not have real roots. For the second condition, these roots are negative, and the corresponding ξ are not real.

Lemma 2. If c ̸ = 0, then the determinant of the matrix M (ξ) is different from 0 for all real ξ if and only if both of the following conditions are satisfied:

where

). The proof of this lemma is straightforward. We note that the real part of det M (ξ) can vanish with its imaginary part different from 0. Then the conditions of the lemma are satisfied and, as it follows from the next theorem, the operator can satisfy the Fredholm property but a part of its essential spectrum can be in the right-half plane resulting in the instability of solutions. We can now prove the Fredholm property of the operator L.

Proof of Theorem 3.2. It follows from the conditions of the theorem, Proposition 3.1 and Lemmas 1, 2 that the operator L is normally solvable with a finite dimensional kernel. It remains to prove that the codimension of its image is finite. Consider an operator L τ such that its dependence on the parameters τ ∈ [0, 1] is continuous in the operator norm, it is normally solvable with a finite dimensional kernel for all τ, L 0 = L and the operator L 1 has a finite codimension of the image. Since the index of the operator L τ does not change during such deformation, then the operator L 0 also has a finite codimension of the image. Thus, in order to prove the theorem, we need to construct a continuous deformation of the operator L to some model operator L 1 in such a way that conditions of Lemma 1 or 2 are satisfied, and the operator L 1 satisfies the Fredholm property. We will reduce the operator L to the operator with zero cross diffusion coefficients (D 2 = D 4 = 0). In this case, the Fredholm property of such operators is known [START_REF] Volpert | Elliptic Partial Differential Equations[END_REF].

Let us begin with the case where c = 0 and condition (A.1) is satisfied. Consider the operator L τ obtained from the operator L if we replace D

Here k is a positive number which will be chosen below. Then L 0 = L and the operator L 1 has zero cross diffusion coefficients. Moreover,


domain: math spi nlin

text:
Reza Abiri

Ian Michael Zeller

Xiaopeng Zhao

Parametric Analysis of Cancer Dynamics: An Evaluation of Environmental Contributing Factors

Mathematical models can be a useful tool to better understand the interaction between physiological parameters that lead to the cancer state from either within and outside the body. One of the most influential effects on the development and continuation of cancer is the effect due to the environment. Ghaffari and Khazaee (2012) developed a time invariant model to demonstrate the onset of a cancerous state through varying environmental factors. The theoretical study followed two identical twin brothers to eliminate genetic variables, and modeled the development of cancer in term of cellular growth rates (a1, a2), cellular carrying capacity (k1, k2), immune system kill rate (d1, d2), external stressors (c), and competition between cell types (b1, b2, g). While the model was successful at showing convergence to either the cancerous or healthy state, a number of improvements are proposed to improve the accuracy and applications of the model. Among these are the introduction of time variation, optimization of input parameters and the addition of a treatment term (D) which fundamentally alters the convergent states of the model. Figure 1 shows a simple schematic of physical meanings of governing equations which are used to describe the healthy and cancerous cell concentrations in an organ based on the subsequent changes. Also, Figure 2 depicts a phase plane of the cancer cells concentration with respect to healthy cells concentration for a healthy person and all the trajectories converge to healthy stable point. Using phase plane analysis, we found that the healthy cell carrying capacity, healthy cell kill rate by the immune system and the healthy vs cancer cell competition were the driving factors for the convergence of the model. Changing these parameters indicated a time frame about which the onset of cancer may occur and that time frame tends to vary between individuals. This element of the model predicted earlier or later cancer onset effectively creating a susceptibility element to the model on an individual basis. For this particular case each parameter was tested for sensitivity while all other parameters were held at their respective values for the birth state. The sensitivity results are shown below in Table 3. As it can be seen that variable (K1) confirms direct relationship in increasing healthy cells and some variables show inverse and the other no changes in this relationship.

Methods Conclusions

The model also showed a successful treatment parameter that promoted convergence to the healthy state at late stages of the model when cancer was more prevalent. Figure 3 Figure 4 shows some graphs that confirm the susceptibility of human to be infected by cancer in different environments.

Introduction Results

Table 3. Dynamic sensitivities have been found for all of the parameters We adopted and extended the model of Ghaffari and Khazaee (2012) to represent an organ model with two stable states (ߝ=0.00001):

.ݍܧ 1: ݔሶ )ݐ( = ܽ ଵ ݔ 1 - ݔ ݇ ଵ -݀ ଵ + ܿ ݔ -ܾ ଵ ݕݔ -ݔܦߝ .ݍܧ 2: ݕሶ )ݐ( = ܽ ଶ ݕ 1 - ݕ ݇ ଶ -݀ ଶ ݕ + ݔܿ -ܾ ଶ ݔ + ݃ ݕ -ݕܦ
where, the healthy and cancerous cell concentrations were denoted by x and y in equations 1 &2, respectively. The model was parameterized as a function of time using a curve fitting method with second order polynomials. The models were tested for sensitivity to determine which parameters have the most substantial effect of the convergence state of the model. This was done by Mathematica software and by holding the initial conditions for the model constant as each parameter was varied and the effect of the convergence state was evaluated. This was indicative of an earlier or later onset of the disease and effectively provides a metric for the influence of the model input parameters. Additionally, a scaled term (D) was added to the model in order to act as a treatment and promote convergence to the healthy state.

Table 1 shows a description for all parameters used in equations 1 & 2. Meanwhile, Table 2 taken from Ghaffari and Khazaee (2012) shows the changes for different parameters during four defined stages of life for a twin brothers. It has been shown that environmental factors can cause a steady state convergence to a cancer state for one of twin brothers (A) living in a polluted environment whereas the other brother (B) who is not living in a polluted environment will show a state convergence to a healthy state. We tested our model by considering the total range of changes in values of these parameters during these four stages. Table 1. The physical meaning of all of the parameters in the governing equations

Bibliography 1 .

 1 Ghaffari, A. and M. Khazaee, Cancer dynamics for identical twin brothers. Theoretical Biology and Medical Modelling, 2012. 9(1): p. 4. 2. Freedman, H. and S. Pinho, Stability criteria for the cure state in a cancer model with radiation treatment. Nonlinear Analysis: Real World Applications, 2009. 10(5): p. 2709-2715. 3. Allavena, P., et al., The inflammatory micro-environment in tumor progression: the role of tumor-associated macrophages. Critical reviews in oncology/hematology, 2008. 66(1): p. 1-9.

  illustrates these convergences.This theoretical model successfully replicated the results from the literature and added additional complexity to present results as a function of time. This effectively changes the convergence state of the model for individual subjects and creates an element of susceptibility for each individual. Treatment was also successfully modeled to promote healthy state convergence. Future work aims to add more complexity to the model and provide empirical physical significance to the input parameters.

Figure 3 .

 3 Figure 3. The addition of treatment changes the convergence from the cancerous to the healthy state at old age. In other words, the stability from cancerous (left) is shifted to a healthy stability (right) by increasing the degree of treatment (D) from left to right. The figure above shows three different degrees of treatment. The first (left) shows D=0, the second (middle) shows D=4 and the third (right) shows D=20 as D increased the stability shifts from a cancerous to a healthy convergence.

Figure 2 .

 2 Figure2. Phase plane and healthy stable condition for a healthy person.

Figure 4 .

 4 Figure 4. The number of cancer cases & deaths reported by WHO. It not only shows different number of cases in different environment, but also confirms that the locations have more occurrence of cases do not have the more mortality. It can verify the positive results of treatment in some locations (source: josefshomperlenblog.org/).

Figure 1 .

 1 Figure 1. The effective parameters in changing of healthy/cancer cell concentration

Table 2 .

 2 The values for parameters for four defined stages of life (copied from [1]).