Theory of possible effects of the Allee phenomenon on refugia of the Hantavirus epidemic
aa r X i v : . [ n li n . AO ] A ug Theory of possible effects of the Allee phenomenon onrefugia of the Hantavirus epidemic
Niraj Kumar, M. N. Kuperman,
1, 2 and V. M. Kenkre Consortium of the Americas for Interdisciplinary Science and Department of Physics and Astronomy,University of New Mexico, Albuquerque, New Mexico 87131, USA Centro At´omico Bariloche and Instituto Balseiro, 8400 S. C. de Bariloche, ArgentinaConsejo Nacional para las Investigaciones Cient´ıficas y T´ecnicas, Argentina
We investigate possible effects of high order nonlinearities on the shapes of infection refugia ofthe Hantavirus epidemic. We replace Fisher-like equations that have been recently used to describeHantavirus spread in mouse populations by generalizations capable of describing Allee effects thatare a consequence of the high order nonlinearities. We analyze the equations to calculate steadystate solutions. We study the stability of those solutions under physical conditions and compareto the earlier Fisher-like case. We consider spatial modulation of the environment and find thatunexpected results appear, including a bifurcation that has not been studied before.
I. INTRODUCTION
Techniques of nonlinear physics and mathematics arefinding much application in recent times in biological andecological systems, which, in turn, are enriching non-linear science. One example is the spread of epidemics[1, 2, 3, 4, 5], in particular the Hantavirus [6] that hasbeen recently modeled [7, 8, 9, 10] with the help of theFisher equation [11] with internal states representing in-fection or its absence in the mouse population. Our pur-pose in this paper is to study the consequences of Alleemechanisms (to be described below) in the dynamics ofmice that are the carriers of the Hantavirus. Withoutthese effects, the (Fisher-like) equations are [7, 8], ∂M S dT = bM − cM S − M S MK − aM S M I + D ∂ M S ∂X ,∂M I dT = − cM I − M I MK + aM S M I + D ∂ M I ∂X ,∂M∂T = bM − cM − M K + D ∂ M∂X , (1)where we display a 1 − d system for simplicity. Here, M S and M I are, respectively, densities of susceptibleand infected mice and M = M S + M I is the totalmouse density, a controls transmission of infectionon encounter, b and c are rates of birth and deathrespectively and D is the mouse diffusion coefficient.The third of the Eqs.(1) is obtained simply by addingthe first two and is of the standard Fisher form [1].The first two possess the characteristic features of theHantavirus [6] that mice are never born infected andthat they are unaffected in any other way (for instance,they do not die faster) when infected. The logisticreaction term in the last of the Eqs.(1) is made up of alinear growth term ( b − c ) M and bilinear depletion term − M /K . Such a description of population dynamics iswidespread [1]. The steady state homogeneous solutionsfor M are [7]: 0 , K ( b − c ). The first of these is unstablewhereas the second, M = K ( b − c ) is stable. The Fisher analysis of the Hantavirus is based on theuse of a logistic reaction term, whereas the investigationwe report in the following assumes a Nagumo termwhich provides an additional zero in the nonlinearityrelative to the logistic case. The physical content behindsuch a term is the Allee effect , in the presence of which,unlike in the logistic case, the zero- M solution is stable.If M is small initially, it is attracted to the vanishingvalue; if large, it is attracted to the nonzero value. Thephysical origin of the Allee effect is the possible increaseof survival fitness as a function of population size forlow values of the latter. Existence of other membersof the species may induce individuals to live longerwhereas low densities may, through loneliness, lead toextinction. There is a great deal of evidence for suchan effect in nature [12, 13] and there have been recentreports [14, 15] of theoretical work addressing the effect.This paper is set out as follows. The Allee effect isdescribed by adding cubic terms to the logistic depen-dence and the model is displayed in section 2, both indimensioned and dimensionless forms. The former is im-portant to understand the connection of the parametersto quantities observed in nature, while the latter facili-tates mathematical analysis. We also define two quanti-ties, α and χ , important to our later development. Theformer is central to the classification of regimes of behav-ior and is directly associated with the Allee phenomenon.We present a linear stability analysis of the steady statesolutions in section 3 and argue the importance of thethreshold α = 0 . χ exceeds a critical value χ c , we carry out numer-ical solutions of the model equations in section 4 for thesimple case of spatially constant α , and in section 5 forthe richer case of spatially modulated α . Conclusions ap-pear in section 6. II. THE MODEL
Our starting set of equations is, instead of (1), ∂M S dT = bM − cM S − M S M K − aM S M I + D ∂ M S ∂X ,∂M I dT = − cM I − M I M K + aM S M I + D ∂ M I ∂X ,∂M∂T = − cM + bM − M K + D ∂ M∂X . (2)The birth term is quadratic in the mouse density andthe environment population term is cubic. This is incontrast to the Fisher case (1) where these terms are lin-ear and quadratic respectively. Needless to say, using aquadratic versus linear birth rate has nothing to do withsexual versus asexual reproduction. Surely, the mice doreproduce sexually. The effective power to be used inthe equations is the result of a variety of factors influenc-ing one another including the probability of encounter ofmates. The appropriate power of the density can be de-termined only phenomenologically, after the fact. Similarconsiderations apply to the quadratic versus cubic envi-ronment term. The important distinguishing feature ofthe Allee case is the possibility of extinction that smallenough populations have to face as a result of the non-linearity, and mathematically speaking, the existence oftwo (rather than one) stable fixed points.The simplest way to analyze the steady states of thesystem is to first look at the solutions of M from thelast of Eq. (2). There are three possible homogeneoussteady state solutions for M : M = 0 , M mx = ( bK + √ b K − cK ) / M mn = ( bK − √ b K − cK ) / M and M mx are stable while M mn is unstable. However, it is important to note that if b < c/k , there is only one real solution correspond-ing to M = 0. It is convenient to reduce Eq. (2) todimensionless form by performing the following substitu-tions: x = X p M mx /KD,t = M mx T /K,m = M/M mx ,m s = M S /M mx ,m i = M I /M mx . (3)We introduce two new quantities that will prove impor-tant in the subsequent analysis: α = M mn M mx = bK − √ b K − cKbK + √ b K − cK , (4) χ = (1 + α ) ab . (5)Using b/M mx = ( α + 1) /K and c/M mx = α/K , we get the following transformed dimensionless equation set. ∂m s ∂t = ( α + 1) m − αm s − m s m − χm s m i + ∂ m s ∂x ,∂m i ∂t = − αm i − m i m + χm s m i + ∂ m i ∂x ,∂m∂t = m (1 − m )( m − α ) + ∂ m∂x . (6) III. LINEAR STABILITY ANALYSIS
Given that the definition of α in Eq. (5) as a mousedensity ratio excludes consideration of values of α thatare not real, the only steady values for m are: 0, α and1, and there is only one steady stable state for the totalpopulation, corresponding to extinction. The solutions m = 0 , m = α is unstable. The crucialconsequence of the Nagumo term (cubic nonlinearity) isthat the solution m = 0 is stable whereas in the Fishercase it is necessarily unstable. Thus, in the presence ofAllee effects, it is possible that the mouse populationmay vanish completely. Since the system in its steadystate adopts the values of the stable solutions, we needto analyze the mouse density when m = 0 or m = 1.The case m = 0 is trivial. When m = 1, we haveonly two valid solutions for ( m s , m i ), (1 ,
0) and (( α +1) /χ, ( χ − − α ) /χ ). We can observe, by performing alinear stability analysis, that when the second solutionadopts negative values, it is unstable while the first oneis stable. When χ = 1 + α , there is a transition, thefirst solution becoming unstable and the second becomingpositive definite and stable. Thus, 1 + α represents thecritical value of χ : χ c = 1 + α. The results can be summarized as follows:( m s , m i ) = (0 , , (1 ,
0) if α > χ − , , (( α + 1) /χ, ( χ − − α ) /χ ) if α ≤ χ − spatial density profiles of mice. Here, we would liketo stress that in the presence of Allee effects, the infectedmouse density vanishes for a < b . In their absence onthe other hand, when the Fisher equation is appropriate[7], the condition involves not only a and b but also c and K . As χ crosses the value χ c , there is a change inthe stable character of the solutions. The system showsan imperfect pitchfork bifurcation : the number of fixedpoints is two for χ < χ c and three for χ > χ c . Thebifurcation is imperfect as the system under study isnot symmetric under reflection, i.e., m s → − m s and m i → − m i . By contrast, the Fisher modeling of theHantavirus epidemic shows a transcritical bifurcation [7].Our study of this system shows that α = 0 . α . To see this,recast Eq.(6) in the Ginzburg-Landau form, ∂m∂t = − δ F δm , where the functional representing the free energy densityis given by F ( x, t ) = Z dx " (cid:18) ∂m∂x (cid:19) + F ( m ) (8) F ( m ) = − Z m f ( m ′ ) dm ′ . (9)The dynamics is determined by F in that the systemevolves such that d F dt ≤
0. The expression for F containstwo terms. The first term involving the derivative tries tominimize F by minimizing the fluctuations in the density.The minima of second term, i.e., F ( m ) are governed bythe reaction term, f ( m ) = m (1 − m )( m − α ) and so wehave, F ( m ) = m / αm / − (1 + α ) m / . The significance of the threshold value α = 0 . F (0) = F (1), i.e., solutions corresponding to both m = 0and m = 1 are equally stable. When α < . , the solution m = 1 is relatively more stable than the solution m = 0.The situation is reversed for α > . IV. HOMOGENEOUS ENVIRONMENT(CONSTANT α ) Because infected mouse density is always zero irrespec-tive of values of α in the subcritical regime χ < χ c , wefocus our attention below only on steady state mousedensity profiles for the other, more interesting, supercrit-ical regime, χ > χ c . In this case, there is a possibility ofgetting nonzero m i and the steady state depends on theparameter α as well. In order to find the spatial mousepattern, we solve Eqs.(6) numerically with a given ini-tial spatial distribution of susceptible and infected mousedensity.We consider a bounded domain where χ = 2 every-where ( supercritical case since α <
1) and use reflectiveboundary conditions. We analyze various values of α andspatially varying initial densities. Except for α = 0 . α < . m s ( x ) = m i ( x ) = A | cos(2 πωx ) | , a periodically modulated initial conditionfor the initial population, we observe that the system evolves to a steady state characterized by a homogeneoussolution, with m = 1 and a nonzero value for the infectedmouse density, m i = ( χ − − α ) /χ , in the entire spatialdomain. In other words, the system evolves towards thestable state corresponding to m = 1, and not towardsthe other stable state for m = 0. Such a possibility ofchoice between attractors is a feature associated with theAllee effect. It does not arise in the earlier description ofHantavirus spread [7] via Fisher-like equations becausethere there is, in that description, only one attractor. If,on the contrary, α > . , the total mouse density m van-ishes completely in the whole spatial domain. This is, infact, an Allee effect. The observed effect is in contrastto the result obtained using Fisher-like equation wherepopulation can never go to zero because that would cor-respond to an unstable state.The case α = 0 . V. INHOMOGENEOUS ENVIRONMENT(SPATIALLY VARYING α ) Noteworthy effects appear when Allee effects combinewith a spatially varying α . To study them, we maintainthe relationship α ≤ χ − α , viz., the birth parameter b , the death rate c , andthe environment parameter K (see Eqs. (5)), are spa-tially varying and consequently introduce a correspond-ing modulation into α .To understand transparently the physical meaning of apostulated modulation of α as we take below, it is helpfulto consider, only for illustration purposes, the case whenthe quantity 4 c/b K is small, and to expand the α ex-pression in Eq. (5) in powers of that quantity. Then weget α ≈ cb K . (10)Any inhomogeneity in the environment, represented bya spatial variation of c , b or k would be reflected in aspatial modulation of α . A larger death rate, a smallerenvironmental parameter or a smaller birth rate wouldcause a decrease in α. We have analyzed the effect of such a spatial varia-tion in two forms: first by considering a sinusoidal mod-ulation with a characteristic wave number and then byconsidering a less regular behavior of α . First we take α ( x ) = 0 . B sin(2 πωx ), with B < . χ > χ c and avoid α <
0. Unexpected mm i FIG. 1: Mean values of the total (solid) and infected (dashed)population as a function of ω , the wavenumber that describesthe modulation of α ( x ) .behavior emerges. The system undergoes a transition asthe wave number ω of the modulation of α crosses a givencritical value, ω c ≈
10. When ω < ω c , the modulation of α induces a modulation in the population. This occursaround 0 . χ − − α χ for the infected population. But as ω exceeds ω c , despitethe modulation of α , the total population of mice adoptsa homogeneous profile with two different values. Thereare two branches, each one corresponding to the steadysolutions shown in Eq.(7) for ω > ω c . Together withthe absence of modulation in the total population we ob-serve that the infected population survives throughoutthe whole domain. While the total population is homo-geneous, both the infected and susceptible population areoscillatory in space and take values different from thosepreviously calculated.Fig 1 displays graphically the results discussed above.The mean value of the total and infected population isplotted as a function ω , the wavenumber that describesthe modulation of α ( x ). A bifurcation is seen at ω c ≈ α ( x ) > .
5. This can be observed in the two plots of Fig.2, where the spatial profile of both the total and infectedpopulations are shown for values of ω above and belowthe value ω c . Next we provide an example of a more realistic behav-ior of the environment that allows us to show that thiseffect survives even when the modulation is not regular,provided regions exist with oscillations of wavenumberover the critical value. In Fig 3 we plot a situation pre-senting an irregular profile for α (also displayed in thefigure)and the solutions for both branches. mm i a X b FIG. 2: Spatial profile of the total (solid) and infected(dashed) mouse population for (a) ω < ω c and (b) ω > ω c (only the non-null branch is shown). Space is plotted on thex-axis in arbitrary units. VI. CONCLUSIONS
In summary, we have studied what consequences thepossible presence of an Allee phenomenon might havein a population of mice carrying the Hantavirus epi-demic. The Allee features have been incorporated thougha Nagumo, i.e. cubic, term in the nonlinearity. We haveanalytically solved and examined the fixed points of theproblem without diffusion, and numerically investigatedthe steady state profiles in the presence of diffusion. Ourfindings are that, first, as expected, in the presence ofthe Allee phenomenon, the mouse population can van-ish completely. This is not possible when the underlyingequation is of the type used earlier [7]. We also observea dependence from initial conditions that are not presentin the earlier analysis. These effects stem from the exis-tence of two, rather than one, stable solutions. Formallystated, under the effect of the Allee phenomenon, the sys-tem exhibits an imperfect pitchfork bifurcation instead ofthe transcritical bifurcation observed in the earlier Fishercase. The most relevant result is the one obtained by in-troducing environmental spatial inhomogeneities. The mm i a x b FIG. 3: Spatial behavior of the total and infected populationas α (also displayed) varies in space. Each plot correspondsto a different branch. Depicted are α (solid line), the totalmouse population m (dashed line) and the infected mousepopulation m i (dotted line). Space is plotted on the x-axis inarbitrary units. observed effect is better observed in the most abstractcase, when the spatial modulation is sinusoidal. We havefound that there is a bifurcation in the behavior of thesystem. The nature of this bifurcation is more evidentwhen calculating the mean value of the population den-sities. We have shown the existence of a critical valueof the spatial modulation wavenumber ω at which thebehavior of the systems completely changes, displayingbistable behavior that depends on the initial conditions.Later, by taking an hypothetic general situation, we haveshown how this effect operates.For the sake of completeness, we describe briefly theanalysis for the physically (observationally) unimpor-tant but mathematically interesting case when α exactly equals 0.5. The system behaves in a rather complexway. Again, we start with initial conditions of the form m s ( x ) = m i ( x ) = A | cos(2 πωx ) | . The interesting featureis that now both steady states are equally stable. Thebasin of attraction of each state is such that the ampli-tude of the initial modulated condition, plays a relevantrole. There are three different regimes characterized by mm i b X c FIG. 4: The evolution of the initial density profile (a) for m (dashed line) and m i (solid line) towards the steady state(c) for α = 0 . a = 0 .
45. Transient patterns are shownin (b). Here x-axis is the space and y-axis corresponds todensity. Units are arbitrary.
1) the system reaching the homogeneous state m = 1 and m i = ( χ − − α ) /χ in the whole spatial domain, 2) thetotal mouse density m vanishing completely in the wholespatial domain, 3) a steady state being characterized byspatial periodic patterns with m oscillating periodicallybetween 0 and 1. We also observe a similar kind of os-cillation for infected mouse density between m i = 0 and m i = ( χ − − α ) /χ . The two first cases are like thoseencountered already in our analysis above for α > . α = 0 . A , the ampli-tude of the initial condition, decreases the system goesfrom case 1 to case 2 but with an intermediate regimecorresponding to case 3. For α = 0 .
5, we note the exis-tence of two critical values of initial conditions at whichwe observe a transition from oscillating periodic struc-tures to homogeneous structures. The first critical valuecorresponds to a transition from oscillating structure toa homogeneous pattern with m = 0. The second criticalvalue corresponds to transition from the oscillating pat-tern to the homogeneous one with m = 1. Thus, in thethe presence of Allee effects, the mouse density dependsnot only on system parameters but also on the initial dis-tribution , a dependence that is quite impossible in Fisherequation treatments [7].Whether rodents of the kind we are describing do or donot exhibit these various Allee effects we have described isa question for the field biologist to pursue. We hope thatthe interesting consequences that we have predicted the-oretically in this paper will stimulate observational workin this direction. Work related to the present study (butquite different in spirit as well as detail) may be found inrecent papers by Kenkre and Kuperman [14] on bacteriain a Petri dish and more recently in the extensive studiesof Clerc et al. [15] who report a number of insights intoAllee effects on pattern formation. Some of the featuresof the present work that distinguish it from those oth-ers is the existence of infected and susceptible subclassesin the mouse population, from the mathematical pointof view the corresponding bifurcations that occur in thesystem, and from the physical point of view the relevance of the study to the spread of epidemics and formation ofrefugia. In a future publication we will stress a morerealistic 2-d depiction of the refugia.This work was supported in part by the NSF undergrant no. INT-0336343 and by NSF/NIH Ecology of In-fectious Diseases under grant no. EF-0326757. [1] J. D. Murray, Mathematical Biology , 2nd edn, New York,Springer(1993).[2] S. P. Petrovskii and Bai-Lian Li,
Exactly Solvable Mod-els of Biological Invasion , Chapman & Hall, CRC Press,2006.[3] Robert Cantrell and Chris Cosner,
Spatial Ecology viaReaction-Diffusion Equations , Wiley series in Mathemat-ical and Computational Biology, 2003.[4]
Frontiers in Mathermatical Biology. Lecture Notes inBiomathematics
Elements of Mathematical Ecology (CambridgeUniversity Press, Cambridge, U K, 2003).[6] T. L.Yates, J. N. Mills, C. A. Parmenter, T. G. Ksiazek,R. R. Parmenter, J. R. Vande Castle, C. H. Calisher, S.T. Nichol, K. D. Abbott, J. C. Young, M. L. Morrison, B.J. Beaty, J. L. Dunnum, R. J. Baker, J. Salazar-Bravo,and C. J. Peters. Bioscience , 989-998 (2002).[7] G. Abramson and V. M. Kenkre, Phys. Rev. E , 011912(2002).[8] V. M. Kenkre in Patterns, Noise and Interplay of Nonlin-earity and Complexity:Proceedings of the PASI on Mod-ern Challenges in Statistical Mechanics , edited by V. M.Kenkre and K. Lindenberg (AIP New York, 2003).[9] V. M. Kenkre, Physica A , 121-126 (2005).[10] V. M. Kenkre, L. Giuggioli, G. Abramson, and G. Camelo- Neto, Eur. Phys. J. B , 461-470, (2007).[11] R. A. Fisher, Ann. Eugen, London , 355-369(1937).[12] W. C. Allee, The Social Life of Animals ( Beacon press,Boston, 1938).[13] F. Courchamp, T. Clutton-Brock, and B. Grenfell,Trends Ecol. Evol. , 405 (1999); P. A. Stephens andW. J. Sunderland, ibid. , 401 (1999); C. W. Fowlerand J. D. Baker, rep. Int. Whal. Comm. , 545 (1999).[14] V. M. Kenkre and M. N. Kuperman, Phys. Rev. E ,051921(2003).[15] M. G. Clerc, D. Escaff and V. M. Kenkre, Phys. Rev. E72