Spatial Heterogeneity in Epidemiological Models




(1)
Department of Mathematics, University of Florida, Gainesville, FL, USA

 




15.1 Introduction


Classical epidemic models assume that the entire population lives in one area and is well mixed. That assumption is not necessarily satisfied. For instance, the population may be living on isolated islands or in different countries while traveling from one location to another. This spatial heterogeneity affects the transmission of the disease. To understand the precise impact of spatial heterogeneity on the dynamics of a disease, we have to build models that account for that heterogeneity. Spatially explicit models are also more effective in evaluating control strategies, particularly those applied to movement of individuals.

Multiple modeling approaches have been used to account for space and movement of individuals in epidemic models. Space can be incorporated as a variable in discrete or continuous form. One of the most popular discrete-space modeling approaches is the metapopulation approach [16, 52, 102]. We will introduce the metapopulation approach in the next section. Other discrete space modeling approaches include epidemic spatial networks [20, 21, 43, 123], cellular automata [101, 147], and lattice epidemic models [137, 141]. Continuous space models include diffusion epidemic models [34, 82, 156], integrodifferential equation epidemic models [116, 136] and integrodifference equation epidemic models [8].

In this chapter, we will focus on metapopulation models and diffusion models as two types of modeling techniques incorporating movement of individuals. Each of these model types requires a different mathematical tool. Metapopulation models typically represent a large system of ordinary differential equations. Diffusion epidemic models, on the other hand, are small systems of partial differential equations.


15.2 Metapopulation Modeling of Epidemic Spread


The concept of metapopulation does not originate in epidemiology, but in ecology. A metapopulation [95] is a group of populations of the same species that leave in spatially isolated areas but interact on some level. Metapopulations occur when different populations live in fragmented habitats but are connected through migration. The isolated areas that are occupied by each population are called patches. In epidemiology, patches may be cities, countries, islands, or other geographically autonomous regions. A necessary requirement of a metapopulation is that patches be connected through migration. Migration is defined as physical movement of individuals from one area to another. Movement can be short-term or long-term. In short-term movement, individuals visit another location for a period of time and return to the home patch. Even though the movement is short-term, it still allows an infected individual to transmit the pathogen to a susceptible individual of a different patch, thus spreading the disease to other locations. Long-term migration arises when individuals move to another location and settle there. The two types of movement are modeled differently. Short-term movement has been called Lagrangian movement, and the corresponding models are called Lagrangian metapopulation models. Long-term movement has been called Eulerian movement, and the corresponding models are called Eulerian metapopulation models [44].

Metapopulation epidemic models consist of n patches. The population of each patch is assumed to be homogeneously mixing. It is divided into the typical epidemiological classes of susceptibles, infectives, and other classes. The sizes of each of these classes are different on different patches. Individuals of some or all of the classes travel between the patches, which leads to the movement of the disease (Fig. 15.1).

A304573_1_En_15_Fig1_HTML.gif


Fig. 15.1
Schematic representation of a metapopulation epidemic model

Historically, early metapopulation epidemic models incorporated short-term movement between the patches. Correspondingly, we begin by describing the Lagrangian modeling framework in the next subsection. Later, Eulerian models were developed, in which the movement is explicit and occurs at certain rates. We will consider an Eulerian modeling framework in the second subsection of this section.


15.2.1 Lagrangian Movement Epidemic Models


We begin by assuming that the total population occupies n regions or patches. The population of patch i is denoted by N i , and the total population size is $$N = N_{1} +\ldots +N_{n}$$ . The population of the ith region is infected by a pathogen and consists of S i susceptible individuals, I i infected individuals, and R i recovered/immune individuals. We have $$N_{i} = S_{i} + I_{i} + R_{i}$$ . In terms of movement, we assume that the members of each region make short visits to at least some of the other regions and return to their home patch. Furthermore, we shall make the simplifying assumption that all members of each region spend the same amount of time visiting other regions, but that time depends on the visited region. While commuting to other regions, infected visitors can transmit the disease to susceptibles in the visited region, while susceptible visitors can acquire the disease from infected members of the visited region. Long-term migration of the members of this community will not be incorporated in this model. Furthermore, we assume for simplicity that the population size of each region remains constant, that is, births μ i N i are balanced by deaths, where μ i is the death rate in region i. With these assumptions, the model takes the form


$$\displaystyle\begin{array}{rcl} S_{i}'& =& \mu _{i}N_{i} - S_{i}\sum _{j=1}^{n}\beta _{ ij}I_{j} -\mu _{i}S_{i}, \\ I_{i}'& =& S_{i}\sum _{j=1}^{n}\beta _{ ij}I_{j} - (\mu _{i} +\gamma _{i})I_{i},\qquad i = 1,\ldots,n,{}\end{array}$$

(15.1)
where we have omitted the equations for the recovered. Here, γ i represents recovery in the ith patch, and β ij are the transmission rates of infected individuals from patch j to susceptible individuals in patch i. We note again that since we are assuming constant population size in each region, the birth rate equals the death rate. In addition, if we add the equations for S i , I i , and R i , we obtain the following differential equation for the population size of patch i: N i ′ = 0.

We obtain the disease-free equilibrium: $$\mathcal{E}_{0} = (N_{1},0,N_{2},0,\ldots,N_{n},0)$$ . If we imagine that commuting between patches does not occur, that is, if β ij  = 0 for i ≠ j, then we can define a reproduction number for each patch:


$$\displaystyle{\mathcal{R}_{i} = \frac{\beta _{ii}N_{i}} {\mu _{i} +\gamma _{i}}.}$$
If $$\mathcal{R}_{i} < 1$$ , then the disease will disappear in patch i if it is isolated from the metapopulation. If $$\mathcal{R}_{i} > 1$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_15_Chapter_IEq5.gif”></SPAN>, then the disease will persist in patch <SPAN class=EmphasisTypeItalic>i</SPAN> if it is isolated from the metapopulation.</DIV><br /> <DIV id=FPar1 class=
Definition 15.1.

We call patch i a sink if $$\mathcal{R}_{i} < 1$$ . We call patch i a source if $$\mathcal{R}_{i} > 1$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_15_Chapter_IEq7.gif”></SPAN>.</DIV></DIV><br /> <DIV class=Para>We can apply the next-generation approach to compute the reproduction number of the system, but that value is implicit and does not give any insights as to when the disease persists and when it dies out in the metapopulation. Following the next-generation approach, we have<br /> <DIV id=Equ2 class=Equation><br /> <DIV class=EquationContent><br /> <DIV class=MediaObject><IMG alt=

(15.2)
Hence we define the basic reproduction number as $$\mathcal{R}_{0} =\rho (FV ^{-1})$$ , where


$$\displaystyle{ FV ^{-1} = \left (\begin{array}{cccc} \frac{\beta _{11}N_{1}} {\mu _{1} +\gamma _{1}} & \frac{\beta _{12}N_{1}} {\mu _{2} +\gamma _{2}} & \ldots & \frac{\beta _{1n}N_{1}} {\mu _{n} +\gamma _{n}}\\ & &\vdots & \\ \frac{\beta _{n1}N_{n}} {\mu _{1} +\gamma _{1}} & \frac{\beta _{n2}N_{n}} {\mu _{2} +\gamma _{2}} & \ldots & \frac{\beta _{nn}N_{n}} {\mu _{n} +\gamma _{n}}\\ \end{array} \right ). }$$

(15.3)
To derive a condition for disease persistence or extinction, we follow [163]. We order the variables $$(S_{1},\ldots,S_{n},I_{1},\ldots,I_{n})$$ , and we consider the Jacobian of system (15.1). It takes the form


$$\displaystyle{ J = \left (\begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \end{array} \right ). }$$

(15.4)
The Jacobian is block-diagonal, since the matrix A 21 is the zeroth matrix. Hence, its eigenvalues are precisely the eigenvalues of A 11 and A 22. The matrix A 11 is equal to $$\mathrm{diag}(-\mu _{1},\ldots,-\mu _{n})$$ . The matrix A 22 is given by


$$\displaystyle{ A_{22} = \left (\begin{array}{cccc} \beta _{11}N_{1} - (\mu _{1} +\gamma _{1})& \beta _{12}N_{1} & \ldots & \beta _{1n}N_{1}\\ & & \vdots & \\ \beta _{n1}N_{n} &\beta _{n2}N_{n}&\ldots &\beta _{nn}N_{n} - (\mu _{n} +\gamma _{n})\\ \end{array} \right ). }$$

(15.5)
The eigenvalues of A 11 are clearly negative. Consequently, the stability of the disease-free equilibrium depends on the eigenvalues of A 22. The matrix − A 22 has nonpositive off-diagonal entries, and thus it possesses the Z-pattern. Then A 22 will have eigenvalues with negative real parts if and only if − A 22 is an M-matrix. For − A 22 to be an M-matrix, it is sufficient that a vector v > 0 exist such that − A 22 v > 0. Let $$v = (1,1,\ldots,1)^{T}$$ . Then the condition − A 22 v > 0 is equivalent to the following conditions being satisfied:


$$\displaystyle{ (\mu _{i} +\gamma _{i}) - N_{i}\sum _{j=1}^{n}\beta _{ ij} > 0,\quad i = 1,\ldots,n. }$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_15_Chapter_Equ6.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(15.6)</DIV></DIV>The above condition can also be rewritten in the form<br /> <DIV id=Equ7 class=Equation><br /> <DIV class=EquationContent><br /> <DIV class=MediaObject><IMG alt= N_{i}\sum _{j\neq i}^{n}\beta _{ ij},\quad i = 1,\ldots,n. }$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_15_Chapter_Equ7.gif”>

(15.7)
So if (15.7) holds, then the disease will be eliminated. However, (15.7) also implies that $$\max \{\mathcal{R}_{1},\ldots,\mathcal{R}_{n}\} < 1$$ , so in this case, all patches must be sinks. However, if even one patch is a source, that is, if there exists j 0 such that $$\mathcal{R}_{j_{0}} > 1$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_15_Chapter_IEq13.gif”></SPAN>, then condition (<SPAN class=InternalRef><A href=15.7) is violated, and the disease persists in the populations. One nonobvious conclusion of (15.7) is that the disease may persist even if all patches are sinks, that is, $$\max \{\mathcal{R}_{1},\ldots,\mathcal{R}_{n}\} < 1$$ .


15.2.2 Eulerian Movement Epidemic Models


To formulate the Eulerian movement epidemic model, we assume that the population lives on n patches. In Eulerian movement, it is assumed that individuals move to another patch and settle there, becoming a part of the population of the host patch. Further, we assume that within each patch, the population mixes homogeneously and the distribution of the disease is described by a classical SIR endemic model. We denote the migration rates from patch j to patch i by m ij . We assume m ii  = 0. Furthermore, the migration rates of susceptible and infective may be different, so the m ij S denote the migration rates of susceptibles, and the m ij I denote the migration rates of infectives. A more complex model of this form was considered in [17]:


$$\displaystyle\begin{array}{rcl} S_{i}'& =& \mu _{i}N_{i} -\beta _{i}S_{i}I_{i} -\mu _{i}S_{i} -\sum _{j=1}^{n}m_{ ji}^{S}S_{ i} +\sum _{ j=1}^{n}m_{ ij}^{S}S_{ j}, \\ I_{i}'& =& \beta _{i}S_{i}I_{i} - (\mu _{i} +\gamma _{i})I_{i},-\sum _{j=1}^{n}m_{ ji}^{I}I_{ i} +\sum _{ j=1}^{n}m_{ ij}^{I}I_{ j},\qquad i = 1,\ldots,n.{}\end{array}$$

(15.8)
In the above model, we have once again omitted the equation of the recovered. The total birth/recruitment rate into patch i is μ i N i , the mortality rate in patch i is μ i , the transmission rate in patch i is β i , and the recovery rate in patch i is γ i . The total population size is given by $$N(t) = S_{1}(t) + I_{1}(t) +\ldots +S_{n}(t) + I_{n}(t)$$ , and it satisfies the differential equation


$$\displaystyle{N'(t) = 0.}$$
Hence, the total population size remains constant.

To find the disease-free equilibrium, let I i  = 0 for $$i = 1,\ldots,n$$ . The disease-free equilibrium is given by $$\mathcal{E}_{0} = (S_{1},\ldots,S_{n},0,\ldots,0)$$ . Then the disease-free equilibrium satisfies the following system of equations:


$$\displaystyle{\mu _{i}N_{i} -\mu _{i}S_{i} -\sum _{j=1}^{n}m_{ ji}^{S}S_{ i} +\sum _{ j=1}^{n}m_{ ij}^{S}S_{ j} = 0,\qquad i = 1,\ldots,n.}$$

We rewrite this system in the form


$$\displaystyle{\mu _{i}S_{i} +\sum _{ j=1}^{n}m_{ ji}^{S}S_{ i} -\sum _{j=1}^{n}m_{ ij}^{S}S_{ j} =\mu _{i}N_{i}\qquad i = 1,\ldots,n.}$$
The matrix of this system is


$$\displaystyle{ M = \left (\begin{array}{cccc} \mu _{1} +\sum _{ j=1}^{n}m_{j1}^{S}& - m_{12}^{S} &\ldots & - m_{1n}^{S}\\ & &\vdots & \\ - m_{n1}^{S} & - m_{n2}^{S}&\ldots &\mu _{n} +\sum _{ j=1}^{n}m_{jn}^{S}\\ \end{array} \right ). }$$

(15.9)
The matrix M possesses the Z-pattern. Furthermore, for $$v = (1,\ldots,1)^{T}$$ , M T v > 0. Hence, M is an M-matrix. As a result, it has a nonnegative inverse, and therefore the system has a nonnegative solution. We conclude that there is a unique nonnegative disease-free equilibrium.

To compute the reproduction number, we use the next-generation approach. The infective variables are $$I_{1},\ldots,I_{n}$$ . We compute the matrices F and V. The new infections are given by the incidence terms. Hence $$F =\mathrm{ diag}(\beta _{1}S_{1},\ldots,\beta _{n}S_{n})$$ . All remaining terms give V:


$$\displaystyle{ V = \left (\begin{array}{cccc} (\mu _{1} +\gamma _{1} +\sum _{ j=1}^{n}m_{j1}^{I})& - m_{12}^{I} &\ldots & - m_{1n}^{I}\\ & &\vdots & \\ - m_{n1}^{I} & - m_{n2}^{I}&\ldots &(\mu _{n} +\gamma _{n} +\sum _{ j=1}^{n}m_{jn}^{I})\\ \end{array} \right ). }$$

(15.10)
The reproduction number is then defined as $$\mathcal{R}_{0} =\rho (FV ^{-1})$$ . It follows from the results in [159] that if $$\mathcal{R}_{0} < 1$$ , the disease-free equilibrium is locally asymptotically stable, and if $$\mathcal{R}_{0} > 1$$ ” src=”/wp-content/uploads/2016/11/A304573_1_En_15_Chapter_IEq23.gif”></SPAN>, then the disease-free is unstable.</DIV><br /> <DIV class=Para>We define the reproduction number of each patch as<br /> <DIV id=Equ11 class=Equation><br /> <DIV class=EquationContent><br /> <DIV class=MediaObject><IMG alt=

(15.11)
Exercise 15.3 asks you to show that for n = 2, the reproduction number $$\mathcal{R}_{0}$$ computed by the next-generation approach is bracketed between the two patch reproduction numbers.


15.3 Spatial Models with Diffusion


Reaction–diffusion equations are one of the most important tools for modeling spatial movement. They are continuous in both space and time and represent second-order partial differential equations. Reaction–diffusion equations offer the best modeling of spatial propagation that occurs as a result of diffusion. They consist of a diffusion component that is the second-order derivative of the unknown function(s) and a reaction part that resembles the right-hand side of an ordinary differential equation model.


15.3.1 Derivation of Reaction–Diffusion Equations


There are various approaches for the derivation of reaction–diffusion equations. We will use the approach that captures the movement discretely and then obtains a continuous equation by passing to the limit. To simplify matters, we will imagine that the movement occurs in one dimensional direction, or along a line. The description here follows that in [33]. Denote the probability that an individual is at location x at time t by p(x, t). We discretize both space and time. Space is discretized with a step Δ x, and time is discretized with a step Δ t. The probability p(x, t +Δ t) that an individual will be at location x at time t +Δ t is the sum of the probabilities of being at location x −Δ x at time t and moving to the right with probability α or being in position x +Δ x at time t and moving to the left with probability 1 −α, that is,


$$\displaystyle{p(x,t +\varDelta t) =\alpha p(x -\varDelta x,t) + (1-\alpha )p(x +\varDelta x,t).}$$

Subtracting p(x, t) from both sides and dividing by Δ t, we have


$$\displaystyle\begin{array}{rcl} \frac{p(x,t +\varDelta t) - p(x,t)} {\varDelta t} & =& \frac{1} {2\varDelta t}[p(x -\varDelta x,t) - 2p(x,t) + p(x +\varDelta x,t)] \\ & & +\frac{1/2-\alpha } {\varDelta t} [p(x +\varDelta x,t) - p(x -\varDelta x,t)].{}\end{array}$$

(15.12)
At this point, we have to make an assumption for the correlation between Δ t and Δ x. In particular, we assume that $$D = (\varDelta x)^{2}/(2\varDelta t)$$ and $$v = 2(1/2-\alpha )\varDelta x/\varDelta t$$ are constants. Because of these definitions of D and v, they have units distance2/time and distance/time respectively. The constant D is called the diffusion constant. Then the above expression becomes


$$\displaystyle\begin{array}{rcl} \frac{p(x,t +\varDelta t) - p(x,t)} {\varDelta t} & =& D\left [\frac{p(x -\varDelta x,t) - 2p(x,t) + p(x +\varDelta x,t)} {(\varDelta x)^{2}} \right ] \\ & & +v\left [\frac{p(x +\varDelta x,t) - p(x -\varDelta x,t)} {2\varDelta x} \right ]. {}\end{array}$$

(15.13)
Passing to the limit as Δ t → 0 and Δ x → 0, we obtain the following differential equation:


$$\displaystyle{ \frac{\partial p} {\partial t} = D\frac{\partial ^{2}p} {\partial x^{2}} + v\frac{\partial p} {\partial x}. }$$

(15.14)
This is the basic equation for directed diffusion, in which case the probabilities of moving left or right are not equal. If we assume that the probabilities of moving left or right are equal, that is, $$\alpha = 1/2$$ , then v = 0, and we obtain the following simplified equation:


$$\displaystyle{ \frac{\partial p} {\partial t} = D\frac{\partial ^{2}p} {\partial x^{2}}. }$$

(15.15)
This is the basic equation that models random diffusion. We will work primarily with random diffusion models.

Equation (15.14) models only the diffusion part of a population model. To incorporate population growth and disease transmission, we have to have a reaction part, which comes from a regular ordinary differential equation epidemic model. For instance, consider an ordinary differential equation SI model


$$\displaystyle{ \begin{array}{l} S' = -\beta SI \\ I' =\beta SI.\end{array} }$$

(15.16)
Now we assume that susceptible and infected individuals are spatially distributed, S(x, t) and I(x, t), where x is a one-dimensional spatial variable. We add diffusion to this model to obtain the following epidemic reaction–diffusion model:


$$\displaystyle{ \begin{array}{l} \frac{\partial S} {\partial t} = -\beta S(x,t)I(x,t) + D\frac{\partial ^{2}S} {\partial x^{2}}, \\ \frac{\partial I} {\partial t} =\beta S(x,t)I(x,t) + D\frac{\partial ^{2}I} {\partial x^{2}},\end{array} }$$

(15.17)
where β is still a constant independent of space. We denote by N(x, t) the total population size. Hence, we have $$N(x,t) = S(x,t) + I(x,t)$$ . The total population size satisfies an equation of the form (15.15):

Models of the form (15.14) and (15.17) will have a unique solution if two types of additional conditions are specified for them:

1.

Initial conditions. These specify given spatial distributions at time t = 0. These are similar to the initial conditions specified for ordinary differential equation models.

 

2.

Boundary conditions. These specify the values of the unknown function(s) at the boundaries of the spatial region.

 

Boundary conditions can be of different types and specify different types of problems when coupled with the same system of differential equations. If the domain is infinite, say the whole real axis, then typically a condition governing the growth of the functions at $$\pm \infty $$ are specified. In particular, with problem (15.17) we will have


$$\displaystyle{S(x,t) \rightarrow 0,\qquad I(x,t) \rightarrow 0\quad \mathrm{as}\quad x \rightarrow \pm \infty.}$$

If the spatial domain is finite, say [0, L], then the following types of boundary conditions are often used:



  • Dirichlet Boundary Conditions. For Dirichlet boundary conditions, the unknown function(s) are prescribed given known values at the boundary of the domain. In particular, for model (15.17), we have


    $$\displaystyle{S(0,t) = S_{0}(t)\qquad S(L,t) = S_{L}(t),\qquad I(0,t) = I_{0}(t),\qquad I(L,t) = I_{L}(t),}$$
    where S 0(t), S L (t), I 0(t), I L (t) are given known functions. A commonly used set of Dirichlet boundary conditions prescribes the values at the boundary to be zero:


    $$\displaystyle{S(0,t) = 0,\qquad S(L,t) = 0,\qquad I(0,t) = 0,\qquad I(L,t) = 0.}$$
    In this case, the Dirichlet boundary conditions are called homogeneous Dirichlet boundary conditions. These conditions are often called “absorbing” or “deadly,” because the boundary “absorbs” all individuals who encounter it.


  • Neumann Boundary Conditions. For the Neumann boundary conditions, the fluxes through the boundary of the unknown function(s) are prescribed given known values of the fluxes at the boundary of the domain. In particular, for model (15.17), we have


    $$\displaystyle{\frac{\partial S} {\partial x} (0,t) =\hat{ S}_{0}(t),\!\!\!\!\!\qquad \frac{\partial S} {\partial x} (L,t) =\hat{ S}_{L}(t),\!\!\!\!\!\qquad \frac{\partial I} {\partial x}(0,t) =\hat{ I}_{0}(t),\!\!\!\!\!\qquad \frac{\partial I} {\partial x}(L,t) =\hat{ I}_{L}(t),}$$
    where $$\hat{S}_{0}(t),\hat{S}_{L}(t),\hat{I}_{0}(t),\hat{I}_{L}(t)$$ are given known functions. Commonly used Neumann boundary conditions are those such that the prescribed flux values at the boundary are zero:


    $$\displaystyle{\frac{\partial S} {\partial x} (0,t) = 0,\qquad \frac{\partial S} {\partial x} (L,t) = 0,\qquad \frac{\partial I} {\partial x}(0,t) = 0,\qquad \frac{\partial I} {\partial x}(L,t) = 0.}$$
    In this case, the Neumann boundary conditions are called homogeneous Neumann boundary conditions. These conditions are often called “no-flux” conditions, because the boundary does not allow individuals to pass through it. These are the most commonly used in population models.


  • Robin Boundary Conditions. For the Robin boundary conditions, a linear combination of the values and the fluxes through the boundary of the unknown function(s) are prescribed given known values at the boundary of the domain. In particular, for model (15.17), we have


    $$\displaystyle\begin{array}{rcl} \alpha _{1}\frac{\partial S} {\partial x} (0,t) +\gamma _{1}S(0,t) =\hat{ S}_{0}(t),& \qquad & \alpha _{2}\frac{\partial S} {\partial x} (L,t) +\gamma _{2}S(L,0) =\hat{ S}_{L}(t), \\ \alpha _{3} \frac{\partial I} {\partial x}(0,t) +\gamma _{3}I(0,t) =\hat{ I}_{0}(t),& \qquad & \alpha _{4} \frac{\partial I} {\partial x}(L,t) +\gamma _{4}I(L,t) =\hat{ I}_{L}(t), {}\end{array}$$

    (15.18)
    where $$\hat{S}_{0}(t),\hat{S}_{L}(t),\hat{I}_{0}(t),\hat{I}_{L}(t)$$ are given known functions. Commonly used Neumann boundary conditions are those such that the prescribed flux values at the boundary are zero:


    $$\displaystyle\begin{array}{rcl} \alpha _{1}\frac{\partial S} {\partial x} (0,t) +\gamma _{1}S(0,t) = 0,& \qquad & \alpha _{2}\frac{\partial S} {\partial x} (L,t) +\gamma _{2}S(L,0) = 0, \\ \alpha _{3} \frac{\partial I} {\partial x}(0,t) +\gamma _{3}I(0,t) = 0,& \qquad & \alpha _{4} \frac{\partial I} {\partial x}(L,t) +\gamma _{4}I(L,t) = 0. {}\end{array}$$

    (15.19)
    In this case, the Robin boundary conditions are called homogeneous Robin boundary conditions. These conditions are often called “mixed” conditions.

A combination of a differential equation and Dirichlet boundary conditions is called a Dirichlet problem. Similarly, a combination of a differential equation and Neumann boundary conditions is called a Neumann problem.


15.3.2 Equilibria and Their Local Stability


Reaction–diffusion equations, just like ordinary differential equations, have time-independent solutions, called equilibria. Equilibrium solutions may be constant in space, in which case they are called spatially homogeneous, or they may depend on the spatial variable explicitly, in which case they are called spatially heterogeneous.

To illustrate this concept, we consider model (15.17) with homogeneous Neumann boundary conditions. The total population size N satisfies the following problem:


$$\displaystyle\begin{array}{rcl} \frac{\partial N} {\partial t} & =& D\frac{\partial ^{2}N} {\partial x^{2}}, \\ \frac{\partial N} {\partial x} (0,t)& =& \frac{\partial N} {\partial x} (L,t) = 0.{}\end{array}$$

(15.20)
If the initial conditions are constant in space, this system clearly has the constant N as a solution, where the value of N is given by the initial conditions. Expressing $$S(x,t) = N - I(x,t)$$ and eliminating S, we may reduce the system (15.17) to the single differential equation


$$\displaystyle\begin{array}{rcl} \frac{\partial I} {\partial t} & =& \beta I(x,t)(N - I(x,t)) + D\frac{\partial ^{2}I} {\partial x^{2}}, \\ \frac{\partial I} {\partial x}(0,t)& =& \frac{\partial I} {\partial x}(L,t) = 0. {}\end{array}$$

(15.21)
Equation (15.21) was first proposed by Fisher [62] to model gene spread in a population and has been widely studied since then. Now this equation is referred to as the Fisher–Kolmogorov equation.

Equilibrium solutions of model (15.21) satisfy


$$\displaystyle\begin{array}{rcl} 0& =& \beta I(x)(N - I(x)) + D\frac{\partial ^{2}I} {\partial x^{2}}, \\ \frac{\partial I} {\partial x}(0)& =& \frac{\partial I} {\partial x}(L) = 0. {}\end{array}$$

(15.22)
This system clearly has two spatially homogeneous solutions: I = 0 and I = N. The system (15.21) does not have spatially heterogeneous solutions. To see this, integrate the first equation in system (15.22) to obtain


$$\displaystyle{\frac{\partial I} {\partial x} = - \frac{\beta } {D}\int _{0}^{x}I(y)(N - I(y))dy.}$$
The first boundary condition is clearly satisfied. However, for nonnegative I, smaller than N, $$\frac{\partial I} {\partial x}$$ is a function that is increasing in absolute value and negative. Thus, $$\frac{\partial I} {\partial x}(L) < 0$$ , which contradicts the second boundary condition.

Local stability of these equilibria is determined by the following theorem:


Theorem 15.1.

The equilibrium I = 0 is unstable. The equilibrium I = N is locally asymptotically stable.


Proof.

We begin with the equilibrium I = 0. We linearize around this equilibrium. In particular, we set I(x, t) = y(x, t). We obtain the following linear equation for the perturbation:


$$\displaystyle\begin{array}{rcl} \frac{\partial y} {\partial t} & =& \beta Ny(x,t) + D\frac{\partial ^{2}y} {\partial x^{2}}, \\ \frac{\partial y} {\partial x}(0,t)& =& \frac{\partial y} {\partial x}(L,t) = 0.{}\end{array}$$

(15.23)
A traditional approach to solving Eq. (15.23) is to look for solutions in separable form, that is, we look for a solution of the form
Get Clinical Tree app for offline access