Bistable reaction-diffusion on a network

Bistable reaction-diffusion on a network

J.-G. Caputo, G. Cruz-Pacheco and P. Panayotaros
Laboratoire de Mathematiques, INSA de Rouen,
Av. de l’Universite, 76801 Saint-Etienne du Rouvray, France
Depto. Matemáticas y Mecánica, I.I.M.A.S.-U.N.A.M.,
Apdo. Postal 20–726, 01000 México D.F., México

We study analytically and numerically a bistable reaction-diffusion equation on an arbitrary finite network. We prove that stable fixed points (multi-fronts) exist for any configuration as long as the diffusion is small. We also study fold bifurcations leading to depinning and give a simple depinning criterion. These results are confirmed by using continuation techniques from bifurcation theory and by solving the time dependent problem near the treshold. A qualitative comparison principle is proved and verified for time dependent solutions, and for some related models.

1 Introduction

Discrete reaction-diffusion equations arise in many different fields. For example they can describe the propagation of a nerve impulse in a neuron [1] or the motion of a dislocation [2]. The solutions of these equations are typically fronts connecting two regions of constant value, say 0 and 1. Front pinning and propagation has been studied by many authors for a one dimensional network for a bistable cubic reaction term. An important result obtained by Keener[13] is that when the Laplacian is weak, any arbitrary configuration of 0’s and 1’s leads to a stable static solution. The study was extended by Erneux and Nicolis[7] who explicitely calculated these fronts and gave a pinning criterion. For material science applications and in the presence of an external forcing, Carpio and Bonilla [4] gave pinning conditions and estimated the front speed. For a two dimensional regular lattice, front propagation was studied by Hoffman and Mallet-Paret[8].

The present article considers arbitrary but finite networks, where to our knowledge there are no works. We address specifically this problem and study analytically and numerically static fronts and how they destabilize in an arbitrary finite network (graph). The reaction term we use is the bistable cubic nonlinearity and the diffusion term is the standard graph Laplacian of the network (see e.g. [9]). Throughout the article, we refer to this equation as the Zeldovich model. We introduce and motivate the bistable reaction-diffusion system by considering how an epidemic propagates on a network. To describe how the epidemic front moves on the network, we extend the standard Kermack-McKendrick model (see e.g. [5] for a recent application) to a network and show how it reduces to a discrete Fisher equation. In contrast to the ODE model, the network Kermack-McKendrick model is not commonly used to describe the spread of an epidemic. The Fisher model only describes the propagation phase. The related Zeldovich model we propose is also new but its cubic bistable nonlinearity has a local excitation threshold, which may be a desirable feature for both geographic networks, where the epidemic spreads from one location to another, and agent-based networks, where the disease spreads from one individual to another.

A first result is the existence of static stable fronts for small diffusivity. The argument combines the implicit function theorem (as in the anticontinuous limit used for other lattice problems, see [15]) with small diffusivity asymptotics for the front amplitudes. The proof also uses a suitable definition for the interface between the active and quiescent sites. The statement is analogous to Keener’s result for the integer lattice [13]. We also show that for large diffusivity the only static solutions are spatially homogeneous.

The existence of these fronts depends on the diffusivity, the nonlinearity, and the local excitation threshold parameters of the model. We focus on the dependence of the static fronts on the diffusivity using numerical continuation techniques. The continuation exhibits the fold structure seen in one dimensional studies [7]. For general networks the depinning diffusivity threshold depends on the front configuration, and a static configuration that becomes unstable can be pinned elsewhere. We compute numerically the depinning thresholds for different static solutions and show that they can be predicted accurately by a simple heuristic expression derived for small diffusivity. By solving the time dependant problem, we verify these findings and see how the connectivity of the network affects the propagation of the fronts above the threshold.

We also obtain qualitative comparison results between different solutions of the Zeldovich equation, showing in particular that ”large” fronts involving large regions of 1’s dominate ”small” fronts. Our study also contains comparison results showing that the Fisher equation describes faster front propagations than both the Zeldovich and Kermack-McKendrick equations. These results are also verified numerically. We see also that the Fisher and Kermack-McKendrick fronts propagate at comparable speeds and are much faster that the Zeldovich fronts. Finally we present numerical results for larger local escitation threshold parameters, showing that the static fronts become wider and travel much faster accross the network when they destabilize.
The article is organized as follows. In section 2 we introduce the Zeldovich equation and discuss the other models. Section 3 studies the fixed points of the Zeldovich equation, presenting theretical and numerical continuation results, as well as a depinning criterion. Section 4 describes comparison results between the solutions of the Zeldovich equation, and between solutions of the Zeldovich, Fisher, and Kermack-McKendrick equations Section 5 presents numerical results of the evolution problem; there we validate the pinning thershold for different fronts and compare the dynamics of large and small fronts. We also show that fronts become wider as the nonlinearity treshold increases and we compute the pinning treshold. Conclusions are given in section 6.

2 The Zeldovich model and epidemic propagation

One of the main models to describe the time evolution of the outbreak of an epidemic is the Kermack-McKendrick model[10]


where are respectively the number of people susceptible to be infected, the number of infected and the number of recovered in a total constant population . We have of course

The dynamics of the model is that (resp. ) if (resp. ) . We also can compute the ”final” state of after the outbreak

Roughly speaking, assuming that is near zero, and , the infected population increases, reaches a maximum value and decreases to zero. The main questions are that maximum value of , the time to reach its, the integral of , etc.

We rescale the variables by

This yields the system


We introduce now the possibility of dispersion from city to city with a Laplacian term. The system (4) becomes


where and the third equation is omitted because of the conservation


This model will describe the outbreak of the epidemic, its spreading, and eventual demise as peaks and starts decreasing at each site.

To simplify even more the model and get analytical results we only consider the maximum outbreak by eliminating the term and only considering the equation for . If then verifies

so that goes to a constant which we can assume to be 1. Then from (2) for , we get the Fisher equation


This equation has two homogeneous solutions and the former is unstable. The model does not have a treshold as opposed to the Kermack-McKendrick. To re-introduce this important feature, we modify the nonlinearity into the cubic (Zeldovich) so that we get


For this, there are only two stable homogeneous solutions . As discussed in the introduction, this equation has many physical applications; it is then an important physical model.

If we had a spatially uniform domain the term would be the usual Laplacian. Here we consider an arbitrary graph, for example the network of six major cities in Mexico shown in Fig. 1. Here the nodes correspond to the cities and the links correspond to the main roads connecting these cities.

Figure 1: Graph of the six main cities in Mexico numbered from 1 to 6: Guadalajara, Zacatecas, Queretaro, Pachuca, Mexico City, Puebla. The links represent the main roads connecting these cities.

For this particular example, the term is


Note that the graph Laplacian is a non-negative symmetric matrix [9]. We use this property below. In physical units the parameter is


where is a diffusion coefficient and is a typical distance between cities. The typical time for the diffusion is then


At this time we assumed the same diffusion coefficient (weight) for all the links of the network. If a node is more or less remote from its neighbors than the other nodes, then one could modify the weight accordingly. With this generalization, we would still have a positive symmetric graph Laplacian.

Let be the triangle . We have the following result.

Lemma 2.1

The unit cube is invariant under the evolution of the Zeldovich (10) and Fisher equations in . The product of the triangles is invariant under the Kermack-McKendrick (2) system with in .

The lemma follows from Propositions 4.1, 4.4 in section 4 below (these do not use any of the results of section 3). It is also easy to show that the corresponding vector fields at the point inwards at the boundaries.

3 Fixed points of the Zeldovich model

We want to describe a situation where only some nodes are excited; in the epidemic context, it means that some nodes are infected and the rest are susceptible. Only the Zeldovich model (10) has such stable fixed points; these are generalized static “fronts” where some nodes are close to one and the rest close to zero. Therefore, in this section, we concentrate on the fixed points of the Zeldovich model (10). We will clarify the situation for the Fisher model (9) below and show why it is less interesting. For definiteness, throughout this section, we consider the 6 node graph from Fig. 1; it is clear that the results can be extended to an arbitrary finite graph.

The fixed point equation we solve is




is the graph Laplacian of (11), and , . We will examine how the fixed points depend on the coupling parameter .

For , and every partition of the set of nodes into three subsets , , we have a solution of of the form , if , , if , , if . Clearly, these are the only solutions of . An inspection of the Jacobian reveals that when is empty, these solutions are stable. On the other hand if is nonempty these solutions are unstable. The number of unstable direction is the number of sites in . The solutions where is empty are generalizations of the fronts that exist for the one dimensional case, they are the main subject of interest of the article.

3.1 Homogeneous fixed points

Let us now consider the case . The homogeneous fixed points can be analyzed for arbitrary . For that consider the system linearized around the fixed point


where the Jacobian matrix has elements


When the fixed points are homogeneous, has a very simple form, it can be written

respectively for , where is the identity. The matrix is then for some real constant , and is . To study the stability it is then convenient to use the basis of orthogonal eigenvectors of the symmetric matrix [9]

where the eigenfrequencies verify

We write


Plugging the above expression into (18) we get the evolution of the amplitude


for the fixed point . Clearly it is stable for any . In a similar way we can show that is always stable. The fixed point is always unstable since we have an eigenvalue .

3.2 Non homogeneous fixed points

For the non homogeneous fixed points the analysis is not so simple. Let us first consider the case but small. The implicit value theorem implies that each solution of can be continued uniquely, that is, it belongs to a unique smooth one-parameter family of satisfying , , provided that is sufficiently small, see e.g. [17]. The solution of the local branch passing from has the same stability as , for sufficiently small. This follows from the fact that all the solutions are hyperbolic.

The numerical solutions below were obtained using the minpack implementation of Powell’s hybrid Newton method [16]. We start from , solving (14) using Newton’s method and step in . After some , we continue stepping but use the pseudo-arc as a parameter[14] because we anticipate a fold. The linear stability of a solution is computed readily by examining the eigenvalues of at , , i.e.


We see numerically that all solutions of with satisfy , forall . This is also shown in Corollary 3.3 below. As we increase the value of along a branch of solutions continued from an solution , the linear stability remains unchanged, until some , depending on the branch, where we see a fold. The branch is then continued by decreasing , until we reach a different solution of the problem. After the fold the number of stable and stable eigenvalues changes. We observe that when is stable, the branch changes stability at the fold, and is unstable. For example, setting , we see that the unstable solution is connected to the stable solution by a branch that has a fold at . In Fig. 2 we show the value of the component at different values of of the fixed point. The other components start, and finish at the same values.

A similar behavior was observed for all the examples examined, except the spatially homogeneous solutions with , , or . From relation (16) one can see that these exist for all . Based on our numerical observations we conjecture that all inhomogeneous fixed points of the problem (we exclude the spatially homogeneous solutions) belong to branches undergoing a fold bifurcation at some positive value of , i.e. we have branches with folds, connecting pairs of solutions. This conjecture can be checked numerically by continuing all fixed points. From the theoretical point of view we can also show that non-spatially homogeneous fixed points cannot exist for arbitrarily large . We have

Proposition 3.1

There is an , such that all , , that satisfy are of the form , with , , or .

The proof is given at the end of this section. The dynamical importance of will be discussed further in the next section. The general idea is that for all initial conditions () should go to one of the two fixed points , , , as .

Figure 2: Component v.s. for a branch connecting the fixed points and .

An interesting problem is the computation of . One idea is to continue all branches starting at solutions and find the largest value of a fold. This computation would give a lower estimate of , since we can not at present rule out the possibility of fixed points not belonging to these branches. Also it is of interest to see whether we can have a family of fixed points that are stable for arbitrarily close to , e.g. a continuous branch having a fold with change of stability at . To obtain a first estimation of we have examined numerically all branches starting from stable solutions for a fixed value of . There are such branches (we exclude , with , ). These are solutions with . In all (non-spatially homogeneous) cases these solutions are connected to an unstable solution of the problem, with . For , the largest value of the fold coupling is , and is observed for the branch connecting the solutions and .

Note that the solution has only one neighbor. This is read from the Laplacian (11). It is reasonable to expect that the solutions that are the last to exist have the least neighbors. We see from (11) that all other solutions with have more that two neighbors, and it is observed that the corresponding branches undergo folds at smaller values of . For example the branch starting from , with two neighbors by (11), undergoes a fold at , while the branch starting from , with four neighbors, undergoes a fold at . The notion of neighbors can be extended to () solutions with . In such cases we can look for the number of external connections to the set , i.e. the number of points having distance one from . We see that more sites in generally imply lower in the corresponding branch. For example the branch starting from , where has one external connection, undergoes a fold at . This is lower that the value of the fold value of the branch starting from above, with two neighbors but fewer peaks. Comparing the values of for the branches corresponding to and , we also see that complementary solutions , (with ), i.e. ones with , that are disjoint and whose union is the set of all nodes, generally have corresponding branches with different fold values.

The solutions not considered in the above enumeration are expected to correspond to branches of solutions that are linearly unstable. Thus, even if we find a static solution such that , we expect that for , almost all initial conditions of the time dependant system (10) go to either , , , as .

To better understand how depends on the type of front and node connectivity, we develop a simple argument that assumes is small, and that all sites except one that we call have values , or , see subsection 3.3 below. This is consistent with what we see numerically, namely that the node that is destabilized first has value approximatelly , see e.g. Fig. 2. Other sites have values that are much lower. The argument is as follows. Call the value of the node that will first destabilize. Then the equation at for is

where is the number of neighbors of that are at 1 and is the connectivity of . This yields


From the continuation study of the static solutions we have seen that for

Combining this observation with (22) yields the estimate for


This estimate is reported in Table 1 below, together with the found numerically. For relatively small, e.g. for used here, we see excellent agreement.

3.3 Asymptotics of the fixed points

In what follows we show some general results on the profile of the fixed points of (10) for , and small. We estimate the decay of the fixed point profiles away from the sites where the solution is near unity; we also see that we can obtain small asymptotics for at all sites. For instance, we show that the amplitude of the equilibrium at the site is

where is the distance of site from the analogue of the “interface” of the configuration, see Lemma 3.6. Roughly speaking, the interface or “front” of an configuration, defined more precisely below, consists of the sites where the solution jumps from zero to unity. The small asymptotic gives us information on the decay of the as we move away form the sites that are near unity. For sites with value near unity it also tells us that are further away from the interface have values that are much closer to unity.

Proposition 3.2 can be also used to compare small solutions continued from different , see Corollary 3.5 below.

The proof of Proposition 3.2 is based on small expansions

valid for all sites . The idea is to insert these expression into (16) and examine the coefficients of the series. We first obtain a less precise, intermediate statement, Lemma 3.6, using induction on the distance from the “interface” between ones and zeros of the solutions. Proposition 3.2 uses the same strategy, and Lemma 3.6.

The precise statements use the following definitions and notation.

Let denote the sites adjacent to the site . Let . Let denote the distance between the set of sites , and a node .

Given a nontrivial solution of the equation , denote by , , the sets of indices where , , respectively. Also let . Let be the set of nodes having at least one neighbor . The set plays the role of the “interface” of the configuration.
Then we have:

Proposition 3.2

Let be a nontrivial solution of equation (14), (16) with , and let , denote the unique branch of solutions of , , that continue for . Consider the sets , , , and corresponding to as defined above, with , nonempty. Then for sufficiently small we have that (i) , imply


and (ii) , imply


An immediate consequence is:

Corollary 3.3

Let be a nontrivial solution equation (16) with , and let , denote the unique branch of solutions of (16), , that continue for . Then for and sufficiently small we have , for all sites .

Proof. For sites we have for sufficiently small. For other sites the statement follows form Proposition 3.2.

Remark 3.4

The above asymptotic is appears to be related to the estimate of in 23, and the assumption that the all sites have values , and . Indeed most sites are seen to be from their values at . Note however that the site also has the value (or ) at , and comes near as approaches . The use of the small asymptotic in justifying 23 is not clear.

Another consequence of Proposition 3.2 is that for sufficiently small there exist pairs of static solutions , of the Zeldovich equation satisfying , . The construction is as follows:

Corollary 3.5

Let , , sufficiently small, be continuations of the fixed points , of the Zeldovich equation satisfying


Then for all sufficiently small we have , .

Proof. We consider the three cases , , and . By (ii) . For we have

for small.

For , Proposition 3.2 yields

with , , and

Therefore for small enough.

For , Proposition 3.2 yields

with , , and

Again for small enough.

The proof of Proposition 3.2 uses the following intermediate result.

Lemma 3.6

Let be a nontrivial solution equation with , and let , denote the unique branch of solutions of , , that continue for . Consider the sets , , , and corresponding to as defined above, with , nonempty. Then for sufficiently small we have that (i) , imply


and (ii) , imply


Proof. We use the analytic version of the implicit value theorem, which allows us to write as a convergent power series in , for sufficiently near the origin. Thus we write , for all sites , see e.g. [17]. (Since the network is finite it is sufficient to use the version for sufficiently large.)

We then already have , , and , .

To show (i) let satisfy . We have


since implies , hence .



By (31), (32), and we must then have .

We use induction: suppose that if , then .

Then for satisfying we have


since implies , hence by the inductive hypothesis. On the other hand


By (33), (34), and we must then have , and therefore , as required.

To see (ii) let satisfy , so that all satisfy . Also . Then


On the other hand


By (35), (36), and we must have , and therefore .

For the inductive step, assume that if satisfies , then . Consider then a site satisfying , then


using the fact that implies , hence by the inductive hypothesis. On the other hand


By (37), (38), implies , and therefore , as required.

We now prove Proposition 3.2.

Proof. The starting point is again the expression . To see (i) first consider sites satisfying . Letting be the set of sites , and , we have



for sufficiently small. On the other hand


By (39), (40), , we need .

We proceed inductively, assuming that if satisfies , then , with . Consider then a site satisfying . Let be the set of sites satisfying . Clearly . Also let . By Lemma 3.6, . Then


for small, since , , by the inductive hypothesis. On the other hand


By (41), (42), and , we therefore need .

To see (ii) consider a site , so that . Then , and


Suppose , then , and (43) yield


for sufficiently small. If , implies , so that (43) implies


for sufficiently small. Combining (44), (45) with


we see that to satisfy with , sufficiently small we must have .

For the inductive step, assume that , imply with . Then let , . Let be the set of sites satisfying , let be the set of sites satisfying .

By Lemma 3.6 we have . Then


for sufficiently small, since , by the inductive hypothesis. On the other hand,


By (47), (48) to satisfy we must have , as required.

We now prove Proposition 3.1

Proof. To study large solutions of we will equivalently examine solutions , where



Then , , is equivalent to , with .

Consider a a sequence , satisfying </