DTIC ADA519129: Maximal Sensitive Dependence and the Optimal Path to Epidemic Extinction

0.54 MB·English

by  Defense Technical Information Center

BulletinofMathematicalBiology(2010)
DOI10.1007/s11538-010-9537-0
ORIGINAL ARTICLE
Maximal Sensitive Dependence and the Optimal Path
to Epidemic Extinction
EricForgostona,∗,SimoneBiancob,LeahB.Shawb,IraB.Schwartza
a
NonlinearSystemsDynamicsSection,PlasmaPhysicsDivision,Code6792,USNaval
ResearchLaboratory,Washington,DC20375,USA
b
DepartmentofAppliedScience,TheCollegeofWilliam&Mary,P.O.Box8795,
Williamsburg,VA23187-8795,USA
Received:13January2010/Accepted:11March2010
©SocietyforMathematicalBiology2010
Abstract Extinctionofanepidemicoraspeciesisarareeventthatoccursduetoalarge,
rarestochasticfluctuation.Althoughtheextinctionprocessisdynamicallyunstable,itfol-
lowsanoptimalpaththatmaximizestheprobabilityofextinction.Weshowthattheopti-
malpathisalsodirectlyrelatedtothefinite-timeLyapunovexponentsoftheunderlying
dynamicalsysteminthattheoptimalpathdisplaysmaximumsensitivitytoinitialcondi-
tions.Weconsiderseveralstochasticepidemicmodels,andexaminetheextinctionprocess
inadynamicalsystemsframework.Usingthedynamicsofthefinite-timeLyapunovex-
ponentsasaconstructivetool,wedemonstratethatthedynamicalsystemsviewpointof
extinctionevolvesnaturallytowardtheoptimalpath.
Keywords StochasticdynamicalsystemsandLyapunovexponents·Optimalpathto
extinction
1. Introduction
Control and eradication of infectious diseases are among the most important goals for
improvingpublichealth.Althoughtheglobaleradicationofadisease(e.g.,smallpox)has
beenachieved(BremanandArita,1980),itisdifficulttoaccomplishandremainsapublic
health target for polio, malaria, and many other diseases, including childhood diseases
(Aylwardetal.,2000).Theglobaleradicationofadiseaseisnottheonlytypeofdisease
extinctionprocess.Forexample,adiseasemay“fadeout”orbecomelocallyextinctina
region.Sincetheextinctionislocal,thediseasemaybereintroducedfromotherregions
(Grassly et al., 2005). Additionally, extinction may also occur to individual strains of a
multistraindisease(MinayevandFerguson,2009),suchasinfluenzaordenguefever.It
isworthnotingthattheextinctionprocessisalsoofinterestinthefieldsofevolutionand
∗
Correspondingauthor.
E-mailaddresses:[e-mail](EricForgoston),[e-mail](SimoneBianco),
[e-mail](LeahB.Shaw),[e-mail](IraB.Schwartz).
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it

does not display a currently valid OMB control number.
1. REPORT DATE 3. DATES COVERED
2010 2. REPORT TYPE 00-00-2010 to 00-00-2010
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Maximal Sensitive Dependence and the Optimal Path to Epidemic
5b. GRANT NUMBER
Extinction
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
Nonlinear Systems Dynamics Section,Plasma Physics Division,Code 6792, REPORT NUMBER
US Naval Research Lab,Washington,DC,20375
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT

15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON
a. REPORT b. ABSTRACT c. THIS PAGE Same as 20
unclassified unclassified unclassified Report (SAR)
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
Forgostonetal.
ecology.Asanexample,bio-diversityarisesfromtheinterplaybetweentheintroduction
andextinctionofspecies(Azaeleetal.,2006;BanavarandMaritan,2009).
In order to predict the dynamics of disease outbreak and spread as well as to im-
plement control strategies which promote disease extinction, one must investigate how
modelparametersaffecttheprobabilityofextinction.Forexample,Dykmanetal.(2008)
showedthatdiseasecontrolandextinctiondependonbothepidemiologicalandsociolog-
icalparametersdeterminingepidemicgrowthrate.Additionally,sinceextinctionoccurs
in finite populations, another factor in determining extinction risk is the local commu-
nitysize(Bartlett,1957,1960;KeelingandGrenfell,1997;ConlanandGrenfell,2007).
Furthercomplicationsarisefromthefactthatingeneral,stochasticeffectscauserandom
transitionswithinthediscrete,finitepopulations.Thesestochasticeffectsmaybeintrin-
sictothesystemormayarisefromtheexternalenvironment.Smallpopulationsize,low
contactfrequencyforfrequency-dependenttransmission,competitionforresources,and
evolutionarypressure(deCastroandBolker,2005),aswellasheterogeneityinpopula-
tionsandtransmission(Lloydetal.,2007)mayallbedeterminingfactorsforextinctionto
occur.Otherfactorswhichcanaffecttheriskofextinctionincludethenatureandstrength
ofthenoise(MelbourneandHastings,2008),diseaseoutbreakamplitude(Alonsoetal.,
2006),andseasonalphaseoccurrence(Stoneetal.,2007).
Inlargepopulations,theintensityoftheintrinsicnoiseisoftenquitesmall.However,
itispossiblethatarare,largefluctuationwhichoccurswithfiniteprobabilitycancause
thesystemtoreachtheextinctstate.Sincetheextinctstateisabsorbingduetoeffective
stochastic forces, eventual extinction is guaranteed when there is no source of reintro-
duction(Bartlett,1949;AllenandBurgin,2000;Gardiner,2004).However,becausefade
outsareusuallyrareeventsinlargepopulations,typicaltimescalesforextinctionmaybe
extremelylong.
Birth-death systems (Gardiner, 2004; van Kampen, 2007), which possess intrinsic
noise, form an important class of stochastic processes. These systems have been used
inthefieldofmathematicalepidemiology(Bartlett,1961;AnderssonandBritton,2000;
Choisyetal.,2007).Inthesesystems,theintrinsicnoisearisesfromthediscretenessof
theindividuals(orspecies)andtherandomnessoftheirinteractions.Topredictprobabili-
tiesofeventsoccurringatcertaintimes,adescriptionofthestochasticsystemisprovided
usingthemasterequationformalism.Stochasticmasterequationsarecommonlyusedin
statisticalphysicswhendealingwithchemicalreactionprocesses(Kubo,1963).
Forsystemswithalargenumberofindividuals,avanKampensystem-sizeexpansion
maybeusedtoapproximatethemasterequationbyaFokker–Planckequation(Gardiner,
2004;vanKampen,2007).However,thetechniquefailsindeterminingtheprobabilityof
largefluctuations(Gaveauetal.,1996;ElgartandKamenev,2004).Instead,inthisarticle,
wewillemployaneikonalapproximationtorecasttheproblemintermsofaneffective
classicalHamiltoniansystem(KamenevandMeerson,2008).Withhighprobability,the
intrinsicnoisewillinduceextinctionofthediseaseorspeciesalongaheteroclinictrajec-
toryinthephasespaceoftheclassicalHamiltonianflow.Thistrajectoryisknownasthe
optimalpath(toextinction).
It is highly desirable to locate the optimal path since the extinction rate depends on
theprobability to traverse thispath. Additionally,the effect of a controlstrategy on the
extinctionratecanbedeterminedbyitseffectontheoptimalpath(Dykmanetal.,2008).
Throughtheuseoftheeikonalapproximation,weconsideramechanisticdynamicalsys-
temsproblemratherthantheoriginalstochasticproblem.Unlikeothermethods,thisap-
proach enables one both to estimate extinction lifetimes and to draw conclusions about
MaximalSensitiveDependenceandtheOptimalPathtoEpidemic
thepathtakentoextinction.Thismoredetailedunderstandingofhowextinctionoccurs
mayleadtonewstochasticcontrolstrategies(Dykmanetal.,2008).
Ingeneral,theoptimalpathtoextinctionisanunstabledynamicalobject.Therefore,
many researchers have investigated the scaling of extinction rates with respect to a pa-
rameter near a bifurcation point, where the dynamics are slow (Doering et al., 2005;
KamenevandMeerson,2008;Kamenevetal.,2008;Dykmanetal.,2008).Theanalytical
calculationofextinctionratesfarfromabifurcationpointismuchmoredifficult,andthus
researchersoftenresorttoaveragingovermanystochasticruns(e.g.,ShawandSchwartz,
2010).Thenumericalcomputationoftheoptimalpathtrajectoryhasbeenachievedinthe
pastusingashootingmethod(KamenevandMeerson,2008).However,sincetheproce-
dureisverysensitivetoboundaryconditions,itisdifficulttoimplementwhenanalyzing
pathsfarawayfrombifurcationpointsorwhenanalyzinghigh-dimensionalmodels.
In this article, we develop a novel method for computing the optimal extinction tra-
jectorythatreliesonthecalculationofthedynamicalsystem’sfinite-timeLyapunovex-
ponents (FTLE). The classical Lyapunov exponent provides a quantitative measure of
howinfinitesimallycloseparticlesinadynamicalsystembehaveasymptoticallyastime
t →±∞ (Guckenheimer and Holmes, 1986). Similarly, the FTLE provides a quantita-
tivemeasureofhowmuchnearbyparticlesseparateafteraspecificamountoftimehas
elapsed.
The FTLE fields were used by Pierrehumbert (1991) and Pierrehumbert and Yang
(1993)tocharacterizestructures,includingmixingregionsandtransportbarriers,inthe
atmosphere. Later, these structures were quantified more rigorously by introducing the
ideaofLagrangianCoherentStructures(LCS)(Haller,2000,2001,2002;Shaddenetal.,
2005; Lekien et al., 2007; Branicki and Wiggins, 2010). The definition of a LCS as a
ridgeoftheFTLEfieldwasintroducedbyHaller(2002)andformalizedbyShaddenetal.
(2005).TheFTLEmethodprovidesameasureofhowsensitivelythesystem’sfuturebe-
haviordependsonitscurrentstate,andtheLCS,orridge,isastructurewhichhaslocally
maximal FTLE value. In this article, we will show that the system displays maximum
sensitivityneartheoptimalextinctiontrajectory,whichenablesustodynamicallyevolve
towardtheoptimalescapetrajectoryusingFTLEcalculations.
In this article, we consider three standard epidemic models that contain intrinsic or
externalnoisesourcesandillustratethepowerofourmethodtolocatetheoptimalpath
to extinction. Even though our examples are taken from infectious disease models, the
approach is applicable to any extinction process or escape process. Section 2 discusses
stochastic modeling in the limit of large population size, while Section 3 describes the
theorythatunderliestheLyapunovexponentcomputations.InSection 4,ourmethodis
usedtofindtheoptimalpathtoextinctionforthreeexamples.InSection5,wedemon-
strate that the optimal path to extinction also possesses a local maximum of the FTLE
field,andconclusionsarecontainedinSection6.
2. Stochasticmodelinginthelargepopulationlimit
Tointroducetheideaofconstructinganoptimalpathinstochasticdynamicalsystems,we
considertheproblemofextinctiontakenfrommathematicalepidemiology.Thestochastic
formulation of the problem accounts for the random state transitions which govern the
dynamicsofthesystem.Toquantitativelyaccountfortherandomnessinthesystem,we
Forgostonetal.
willformulateamasterequationwhichdescribesthetimeevolutionoftheprobabilityof
finding the system in a particular state at a certain time (Gardiner, 2004; van Kampen,
2007).
The state variables X∈Rn of the system describe the components of a population,
whiletherandomstatetransitionswhichgovernthedynamicsaredescribedbythetran-
sition rates W(X,r), where r ∈Rn is an increment in the change of X. Consideration
ofthenetchangeinincrementsofthestateofthesystemresultsinthefollowingmaster
equationfortheprobabilitydensityρ(X,t)offindingthesysteminstateXattimet:
∂ρ(X,t) (cid:2)(cid:3) (cid:4)
= W(X−r;r)ρ(X−r,t)−W(X;r)ρ(X,t) . (1)
∂t
r
If the total population size of the system is N, the components of the state variable
canberescaledtobefractionsofthepopulationbyletting x=X/N.Thenthegeneral
solution(Kuboetal.,1973)forthetransitionprobabilityfromx toxinthetimeinterval
0
fromt tot canbewrittenasthefollowingpathintegral:
0
(cid:5) (cid:6) (cid:5) (cid:9)
t (cid:3) (cid:7) (cid:8) (cid:4)
P(x,t|x ,t )= dD(x,p)exp −N ds H x(s),p(s),s −p(s)x˙(s) , (2)
0 0
t0
wheredD(x,p)denotesintegrationoverallpaths.
TheintegralintheexponentofEq.(2)istheaction,andtheHamiltonianH(x,p;t)is
giveningeneralas
(cid:2) (cid:7) (cid:8)
H(x,p;t)= w(x;r) exp(p·r)−1 , (3)
r
where w(x;r) is defined as the transition rate W per individual. The Hamiltonian H
dependsbothonthestateofthesystemxaswellasthemomentump,whichprovidesan
effectiveforceduetostochasticfluctuationsonthesystem.Itshouldbenotedthatinstead
of using the Hamiltonian representation, one could use the Lagrangian representation,
whichresultsinthefollowingalternativesolution:
(cid:5) (cid:6) (cid:5) (cid:9)
t (cid:7) (cid:8)
P(x,t|x ,t )= dD(x)exp N dsL x(s),x˙(s),s , (4)
0 0
t0
withL(x,x˙;t)=−H(x,p;t)+x˙ ·p.
The action reveals much information about the probability evolution of the system,
fromscalingnearbifurcationpointsinnon-Gaussianprocessestoratesofextinctionasa
functionofepidemiologicalparameters(Dykman,1990;Dykmanetal.,2008).Inorder
tomaximizetheprobabilityofextinction,onemustminimizetheaction.Theminimizing
formulation entails finding the solution to the Hamilton–Jacobi equation, which means
that one must solve the 2n-dimensional system of Hamilton’s equations for x and p.
The appropriate boundary conditions of the system are such that a solution starts at a
nonzerostate,suchasanendemicstate,andasymptoticallyapproachesoneormorezero
componentsofthestatevector.Therefore,atrajectorythatisasolutiontothetwo-point
boundaryvalueproblemdeterminesapath,whichinturnyieldstheprobabilityofgoing
MaximalSensitiveDependenceandtheOptimalPathtoEpidemic
from the initial state to the final state. The optimal path to extinction is the path that
minimizestheactionineithertheHamiltonianorLagrangianrepresentation.
To illustrate an application of the theory for a finite population, we consider in de-
tailanexplicitexample,namelythewell-knownproblemofextinctioninaSusceptible-
Infectious-Susceptible(SIS)epidemicmodel.Inthisexample,thepopulationconsistsof
S susceptible individuals and I infectious individuals. The population is driven via the
birthrateμ,whichisalsoequaltothedeathrate.ThetransitionratesforX=(S,I)T are
givenasfollows:
(cid:7) (cid:8) (cid:7) (cid:8)
W X;(1,0) =Nμ, W X;(−1,0) =μX ,
1
(cid:7) (cid:8) (cid:7) (cid:8)
W X;(0,−1) =μX , W X;(1,−1) =γX , (5)
2 2
(cid:7) (cid:8)
W X;(−1,1) =βX X /N,
1 2
whereβ isthemassactioncontactrate,γ istherecoveryrate,andN isnowaparameter
fortheaveragesizeofthepopulation.ForlargeS,I ∝N,fluctuationsofSandI aresmall
onaverage.Ifthesefluctuationsaredisregarded,onearrivesatthefollowingdeterministic
mean-fieldequationsforS andI:
X˙ =Nμ−μX +γX −βX X /N, (6a)
1 1 2 1 2
X˙ =−(μ+γ)X +βX X /N. (6b)
2 2 1 2
Equations(6a)–(6b)arethestandardequationsoftheSISmodelintheabsenceoffluc-
tuations.Fortheparameter R =β/(μ+γ)>1,theyhaveastable,attracting solution
0
X =Nx with x =R−1, and x =1−R−1, which corresponds to endemic dis-
A A 1A 0 2A 0
ease.Inaddition,Eqs.(6a)–(6b)haveanunstablestationarystate(saddlepoint)givenby
X =Nx with x =1 andx =0,whichcorrespondstotheextinct,ordisease-free,
S S 1S 2S
state.
ForN (cid:7)1,thesteadystatedistributionρ(X)hasapeakatthestablestateX with
A
width∝N1/2.ThispeakisformedoveratypicalrelaxationtimegiveninDykmanetal.
(2008)andSchwartzetal.(2009).However,intheprocessofextinction,weareinterested
intheprobabilityofhavingasmallnumberofinfectiousindividuals,whichisdetermined
bythetailofthedistribution.Thedistributiontailcanbeobtainedbyseekingthesolution
ofEqs.(6a)–(6b)intheeikonalform(ElgartandKamenev,2004;Doeringetal.,2005;
Kubo et al., 1973; Wentzell, 1976; Gang, 1987; Dykman et al., 1994; Tretiakov et al.,
2003)givenby
(cid:3) (cid:4)
ρ(X)=exp −NS(x) , x=X/N,
(7)
ρ(X+r)≈ρ(X)exp(−p·r), p=∂S(x)/∂x,
whereS(x)istheaction.
To leading order in N−1, the equation for S(x) has a form of the Hamilton–Jacobi
equation S˙=−H(x,∂ S;t), where S is the effective action, and the effective Hamil-
x
tonian is given by Eq. (3), with w(x;r)=N−1W(X;r) being the transition rates per
individual.TheactionS(x)canbefoundfromclassicaltrajectoriesoftheauxiliarysys-
temwithHamiltonianH thatsatisfythefollowingequations:
x˙ =∂ H(x,p), p˙ =−∂ H(x,p). (8)
p x
Forgostonetal.
Sincethemaximumoftheprobabilityofextinctionisfoundbyminimizingtheaction,
we compute the trajectory satisfying the Hamiltonian system that has as its asymptotic
limits in time the endemic state as t →−∞ and the extinct state as t →+∞. The ac-
tionthenhastheformfromEq.(2)(Wentzell,1976;Gang,1987;Dykmanetal.,1994;
Tretiakovetal.,2003):
(cid:5)
∞
S(x )= p·x˙dt, H(x,p)=0. (9)
S
−∞
InEq.(9),theintegraliscalculatedforaHamiltoniantrajectory(x(t),p(t))T thatstarts
as t →−∞ at x →x ,p→0, and arrives as t →∞ at the state x . This trajectory
A S
describes the most probable sequence of elementary events x →x+r that brings the
systemtox .
S
Several authors have considered how extinction rates scale with respect to a pa-
rameter near bifurcation points (Doering et al., 2005; Kamenev and Meerson, 2008;
Kamenev et al., 2008; Dykman et al., 2008) when the distance to the bifurcation point
issmallandthedynamicsisveryslow.Foranepidemicmodel,thismeansthattherepro-
ductiverateofinfectionisgreaterthanbutveryclosetoone.However,mostrealdiseases
have reproductive rates of infection greater than 1.5 (Anderson and May, 1991), which
translates into faster growthrates fromtheextinct state. Ingeneral, inorder to getana-
lyticscalingresults,onemustcomputetheoptimalpathusingeithertheHamiltonianor
Lagrangianequationsofmotion.However,farfrombifurcationpoints,oneisseldomable
toperformtherequiredanalysisorcomputation.
Additionally, the computation of the optimal path involves the use of a numerical
shootingmethod(KamenevandMeerson,2008).TheHamiltonianorLagrangianrepre-
sentationofann-dimensionaldynamicalsystemliesin2n-dimensionalspace.Therefore,
for even relatively low-dimensional dynamical systems, the use of a shooting method
to find the optimal path to extinction can be quite problematic. In the next section, we
demonstratehowtoevolvenaturallytotheoptimalpathusingadynamicalsystemsap-
proach.
3. Finite-timeLyapunovexponents
Weconsideravelocityfieldv:R2n×I →R2n givenbytheHamiltonianfieldinEq.(8)
whichisdefinedoverthetimeintervalI =[t,t ]⊂Randthefollowingsystemofequa-
i f
tions:
(cid:7) (cid:8)
y˙(t;t,y )=v y(t;t,y ),t , (10a)
i 0 i 0
y(t;t,y )=y , (10b)
i i 0 0
wherey=(x,p)T ∈R2n,y ∈R2n,andt∈I.
0
Such a continuous time dynamical system has quantities, known as Lyapunov expo-
nents,whichareassociatedwiththetrajectoryofthesysteminaninfinitetimelimit.The
Lyapunovexponentsmeasurethegrowthratesofthelinearizeddynamicsaboutthetra-
jectory.Tofindthefinite-timeLyapunovexponents(FTLE),onecomputestheLyapunov
MaximalSensitiveDependenceandtheOptimalPathtoEpidemic
exponents on a restricted finite time interval. For each initial condition, the exponents
provideameasureofitssensitivitytosmallperturbations.Therefore,theFTLEisamea-
sure of the local sensitivity to initial data. For the purpose of completeness, we briefly
recapitulate the derivation of the FTLE. Details regarding the derivation along with the
appropriatesmoothnessassumptionscanbefoundinHaller(2000,2001,2002),Shadden
etal.(2005),Lekienetal.(2007),andBranickiandWiggins(2010).
The integration of Eqs. (10a)–(10b) from the initial time t to the final time t +T
i i
yieldstheflowmapφti+T whichisdefinedasfollows:
ti
φti+T :y (cid:10)→φti+T(y )=y(t +T;t,y ). (11)
ti 0 ti 0 i i 0
ThentheFTLEcanbedefinedas
(cid:10)
1
σ(y,t,T)= ln λ (Δ), (12)
i |T| max
whereλ (Δ)isthemaximumeigenvalueoftherightCauchy–Greendeformationtensor
max
Δ,whichisgivenasfollows:
(cid:11) (cid:12) (cid:11) (cid:12)
dφti+T(y(t)) ∗ dφti+T(y(t))
Δ(y,t,T)= ti ti , (13)
i dy(t) dy(t)
with*denotingtheadjoint.
Foragiveny∈R2n ataninitialtimet,Eq.(12)givesthemaximumfinite-timeLya-
i
punovexponentforsomefiniteintegrationtimeT (forwardorbackward),andprovidesa
measureofthesensitivityofatrajectorytosmallperturbations.
The FTLE field given by σ(y,t,T) can be shown to exhibit “ridges” of local max-
i
imainphasespace.Theridgesofthefieldindicatethelocationofattracting(backward
timeFTLEfield)andrepelling(forwardtimeFTLEfield)structures.Intwo-dimensional
(2D) space, the ridge is a curve which locally maximizes the FTLE field so that trans-
verse to the ridge one finds the FTLE to be a local maximum. What is remarkable is
thattheFTLEridgescorrespondtotheoptimalpathtrajectories,whichisshownheuris-
tically in Section 5. The basic idea is that since the optimal path is inherently unsta-
ble and observed only through many realizations of stochastic experiments, the FTLE
shows that locally, the path is also the most sensitive to initial data. Figure 1 shows a
schematic that demonstrates why the optimal path has a local maximum to sensitivity.
Ifonechoosesaninitialpointoneithersideofthepathneartheendemicstate,thetwo
trajectorieswillseparateexponentiallyintime.Thisisduetothefactthatbothextinctand
endemicstatesareunstable,andtheconnectingtrajectorydefiningthepathisunstableas
well.Anyinitialpointsstartingneartheoptimalpathwillleavetheneighborhoodinshort
time.
4. FindingtheoptimalpathtoextinctionusingFTLE
We now apply our theory of dynamical sensitivity to the problem of locating optimal
paths to extinction for several examples. We consider the case of internal fluctuations,
Forgostonetal.
Fig.1 Schematicshowingthepathfromtheendemicstate(blue)totheextinctstate(red).Theoptimal
pathleavestheendemicpointalonganunstablemanifoldandconnectstotheextinctstatealongastable
manifold.Themagentaandgreendashedlinesrepresenttrajectoriesoneithersideoftheoptimalpath.The
initialstartingdistancebetweentrajectoriesneartheendemicstateexpandsexponentiallyinforwardtime
(shownbythecyanlines).Locally,thisdemonstratesthatthefinite-timeLyapunovmeasureofsensitivity
withrespecttoinitialdataismaximalalongtheoptimalpath.Thisisevidentintheridgesobservedinthe
evolutionoftheexponents.(Colorfigureonline.)
where the noise is not known a priori, as well as the case of external noise, where the
noiseisspecified.Ineachcase,theinteractionofthenoiseandstateofthesystemsbegin
withadescriptionoftheHamiltonian(orLagrangian)todescribetheunstableflow.Then
thecorrespondingequationsofmotionareusedtocomputetheridgescorrespondingto
maximumFTLE,whichinturncorrespondtotheoptimalextinctionpaths(Schwartzet
al.,2010).
4.1. Example1—Extinctioninabranching-annihilationprocess
Foranexampleofasystemwithintrinsicnoisefluctuationswhichhasananalyticalsolu-
tion,weconsiderextinctioninthestochasticbranching-annihilationprocess
A−→λ 2A and 2A−→μ ∅, (14)
whereλandμ>0areconstantreactionrates(ElgartandKamenev,2004;Assafetal.,
2008). Equation (14) is a single species birth-death process and can be thought of as a
simplifiedformoftheVerhulstlogisticmodelforpopulationgrowth(Nåsell,2001).
The stochastic process given by Eq. (14) contains intrinsic noise which arises from
the randomness of the reactions and the fact that the population consists of discrete in-
dividuals. This intrinsic noise, which can generate a rare sequence of events that cause
the system to evolve to the empty state, can be accounted for using a master equation
approach. The probability P (t) to observe, at time t, n individuals is governed by the
n
followingmasterequation:
P˙n= μ2(cid:3)(n+2)(n+1)Pn+2−n(n−1)Pn(cid:4)+λ(cid:3)(n−1)Pn−1−nPn(cid:4). (15)
MaximalSensitiveDependenceandtheOptimalPathtoEpidemic
Using the formalism of Assaf et al. (2008), Eq. (15) is recast as the following exact
evolutionequationforG(ρ,t):
∂G μ(cid:7) (cid:8)∂2G ∂G
= 1−ρ2 +λρ(ρ−1) , (16)
∂t 2 ∂ρ2 ∂ρ
whereGisaprobabilitygeneratingfunctiongivenby
(cid:2)∞
G(ρ,t)= ρnP (t), (17)
n
n=0
andwhereρ isanauxiliaryvariable.
WesubstitutetheeikonalansatzG(ρ,t)=exp[−S(ρ,t)],whereS istheaction,into
Eq. (16) and neglect the higher-order ∂2S/∂ρ2 term. This results in a Hamilton–Jacobi
equationforS(ρ,t).Byintroducingaconjugatecoordinateq=−∂S/∂ρandbyshifting
themomentump=ρ−1,thenonearrivesatthefollowingHamiltonian:
(cid:11) (cid:12)
μ
H(q,p)= λ(1+p)− (2+p)q qp. (18)
2
Hamilton’sequationsarethereforegivenas
∂H (cid:3) (cid:4)
q˙= =q λ(1+2p)−μ(1+p)q , (19a)
∂p
∂H (cid:3) (cid:4)
p˙=− =p μ(2+p)q−λ(1+p) . (19b)
∂q
TheHamiltoniangivenbyEq.(18)hasthreezero-energycurves.Thefirstisthemean-
fieldzero-energylinep=0,whichcontainstwohyperbolicpointsgivenash =(q,p)=
1
(λ/μ,0)andh =(q,p)=(0,0).Thesecondistheextinctionlineq=0,whichcontains
0
anotherhyperbolicpointgivenash =(q,p)=(0,−1).Thethirdzero-energycurveis
2
non-trivialandisgivenasfollows:
2λ(1+p)
q= . (20)
μ(2+p)
ThesegmentofthecurvegivenbyEq.(20)whichliesbetween−1≤p≤0corresponds
toaheteroclinictrajectory.Ast approaches−∞,thetrajectoryexitsthehyperbolicpoint
h alongitsunstablemanifoldandentersthehyperbolicpoint h alongitsstablemani-
1 2
foldast approaches∞.Thisheteroclinictrajectoryistheoptimalpathtoextinctionand
describes the mostprobable (rare) sequence of events which evolves the system from a
quasi-stationarystatetoextinction(Assafetal.,2008).
ToshowthattheFTLEevolvestotheoptimalpath,wecalculatetheFTLEfieldusing
Eqs.(19a)–(19b).Figure2(a)showstheforwardFTLEplotcomputedusingEqs.(19a)–
(19b)forT =6,withλ=2.0andμ=0.5.InFig.2(a),onecanseethattheoptimalpath
toextinctionisgivenbytheridgeassociatedwiththemaximumFTLE.