Speeding up disease extinction with a limited amount of vaccine

0.23 MB·English

by  M. Khasin

Speeding up disease extinction with a limited amount of vaccine
M. Khasin1, M.I. Dykman1 and B. Meerson2
1Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824 USA and
2Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel
Weconsideroptimalvaccinationprotocolwherethevaccineisinshortsupply. Inthiscase,disease
extinctionresultsfromalargeandrarefluctuation. Weshowthattheprobabilityofsuchfluctuation
canbeexponentiallyincreased byvaccination. Forperiodicvaccinationwithfixedaveragerate,the
optimal vaccination protocol is model independent and presents a sequence of short pulses. The
effectofvaccinationcanberesonantlyenhancedifthepulseperiodcoincideswiththecharacteristic
period of the disease dynamics or its multiples. This resonant effect is illustrated using a simple
epidemicmodel. Ifthesystem is periodically modulated,thepulsesmust besynchronized with the
0 modulation, whereas in the case of a wrong phase the vaccination can lead to a negative result.
1 The analysis is based on the theory of fluctuation-induced population extinction in periodically
0 modulated systems that we develop.
2
n PACSnumbers: 87.23.Cc,02.50.Ga,05.40.-a
a
J
I. INTRODUCTION stochastic dynamics of the epidemics. The underlying
1
3 mechanismis the changeof the rate of largefluctuations
Spreading of an infectious disease is a random pro- leading to disease extinction. For a well mixed popu-

] lation, this rate W is usually exponentially small for
E cess. An important source of the randomness is demo- e
a large total population size N ≫ 1, W ∝ exp(−Q)
graphic noise associated with the stochastic character of e
P
with Q ∝ N, [3, 7–14]. We call Q the disease extinc-
the events of infection, recovery, birth, death, etc. In
.
o a large population demographic noise is small on aver- tion barrier. Vaccination changes the value of Q/N. In
i turn, this changes the disease extinction rate exponen-
b age, and the infection spread leads to an endemic state
tially strongly. This effect was previously discussed for
- where a certain fraction of the population stays infected
q vaccination applied at random [13].
for a long time. The disease, however, can disappear as
[
a result of a large rare fluctuation, an unlikely chain of The goal of this paper is to find an optimal way of

2 events where, for example, susceptible individuals hap- administering a limited amount of vaccine which would
v pen to avoid getting infected while infected ones recover maximally increase the disease extinction rate. We find
0 [1–3]. Then, if there is no influx of infected individu- a vaccination protocol that applies for a broad class of
7
als from the outside, the population will be disease-free. epidemic models. Our approach is based on the obser-
1
Suchspontaneousdisappearanceofadiseaseisanexam- vation that, in a large fluctuation that leads to disease
0
. ple of population extinction studied in stochastic popu- extinction, the population is most likely to evolve in a
1
lation dynamics. well-definedway. Itmovesalongthe mostprobablepath
0
Spontaneous extinction is also important for various inthespaceofthedynamicalvariableswhichcharacterize
0
1 physical and chemical reaction systems. This is a conse- different sub-populations [13–15]. Vaccination perturbs
: quence of the underlying similarity of the dynamics that the system as it moves along the optimal path. One can
v
result from short random events, like collisions between think of vaccination as “force” and its effect as “work”
i
X molecules that lead to chemical reactions and interac- done on the system. This work reduces the barrier Q.
r tionsbetweenindividualsthatleadtothediseasespread. Theproblemthenistomaximizethe workforgivencon-
a Extinction can be understood for different types of sys- straints on the vaccine.
tems within the same general formalism, which provides Optimizationofthe effectofvaccinationresemblesan-
a broader scope for the present paper. Moreover, the other problem of optimal control of random systems:
method of optimal controlof extinction that we propose controlling large fluctuations in noise-driven dynamical
can be applied to systems of various types. systems by applying an external field with a given av-
A conventional way of fighting epidemics is via vac- erage intensity [16, 17]. There are, however, important
cination. If there is enough vaccine, the infection can differences, which come from the very nature of the con-
be eradicated “deterministically” by eliminating the en- trol field. Indeed, vaccination only reduces the number
demic state [4]. The amount of available vaccine, how- of susceptible individuals. In other words, as a control
ever,is ofteninsufficient. The vaccinemaybe expensive, field,vaccinationneverchangessign. Then,aswefind,if
or it may be dangerous to store in large amounts, as in the availableamountofvaccineisconstrainedbyagiven
the case of anthrax, or it may be effectively short-lived mean vaccination rate, the optimal vaccination protocol
due to mutations of the infection agent, as for HIV [5] turns out to be model-independent. This applies also to
and influenza [6]. using a limited amount of medications and other situa-
Even where the endemic state may not be eliminated tions where the control field drives the system only in
deterministically,vaccinationcandramaticallyaffectthe one direction.
2
We assume that vaccination is applied periodically in quasi-continuous vector x =X/N, where N is the char-
time. In this case, the optimal protocol is to apply the acteristic total population size, N ≫1. We assume that
vaccine as a sequence of δ-like pulses. The disease ex- the populationdynamics isMarkovian. Itis quite gener-
tinction rate can strongly depend on the period of this allydescribedby the masterequationfor the probability
sequence. Furthermore, the extinction rate can display distribution P(X,t),
exponentially sharp peaks when the pulse sequence pe-
riod is close to the characteristic period of oscillations P˙(X,t) = [W(X−r,r,t)P(X−r,t)
of the system in the absence of fluctuations, or to its Xr
multiples. We illustrate this resonant phenomenon for − W(X,r,t)P(X,t)]. (1)
the Susceptible-Vaccinated-Infected-Recovered (SVIR)
model. Here W(X,r,t) is the rate of an elementary transition
Epidemics often display seasonal modulation [18]. It X→X+r in which the population size changes by r=
is natural to apply a vaccine with period equal to the (r ,r ,...). Examples of such transitions are infection
1 2
modulation period. As we show, there is a qualitative of a susceptible individual as a result of contacting an
difference between the effect of a periodic vaccination in infected individual, recovery of an infected individual or
this case and in the case where seasonal modulation is arrival of a susceptible individual.
absent. For a system with seasonal modulation, an im- We assume that there is no influx of infected individ-
properly applied pulsed vaccination can actually reduce uals into the population. Therefore, there are no tran-
the disease extinction rate and therefore prolong the du- sitions from states where there are no infected to states
ration of the epidemic. The overall effect of the pulsed where infected are present,
vaccinationcriticallydependshereonthephase atwhich
the periodic pulses are applied. W(X,r,t)=0 for XE =0, rE 6=0, (2)
The analysis of periodic vaccination, with and with-
where subscriptE is usedfor the component ofX which
out seasonal variations, necessitates a general formula-
enumerates infected, X ≡I.
tion of the extinction problemin periodically modulated E
Intheneglectofdemographicnoisethepopulationdy-
stochastic populations. We extend the previous results
namicscanbedescribedbythedeterministic(mean-field)
for single-population systems [19, 20] to provide a com-
equation for the population size x¯,
plete extinction theory for multi-population systems in
theeikonalapproximation,andemphasizethedistinction
x¯˙ = rw(x¯,r,t), w(x,r,t)=W(X,r,t)/N. (3)
from the well-understood problem of switching between
metastablestatesinperiodicallymodulatedsystemswith Xr
noise [21].
It immediately follows from Eq. (1) if the width of the
Section II describes the class of epidemic models we
probability distribution P(X,t) is set equal to zero.
consider in this work and develops an eikonal theory
of disease extinction rate in periodically modulated sys-
tems. SectionIIIformulatestheoptimizationproblemfor
1. Stationary systems
vaccination and presents its solution for a time-periodic
vaccination in the limit of a small average vaccination
We start with the case where the transition rates
rate. In Section IV we discuss the vaccination-induced
W(X,r) are independent of time. An endemic state,
reduction of the disease extinction barrier for two types
whereafinitefractionofpopulationisinfectedforalong
of constraints on the vaccination period, a limited life-
time, corresponds to an attracting fixed point x of the
time of the vaccine and a limited vaccine accumulation. A
dynamical system, Eq. (3). We will assume throughout
Section V illustrates, on the example of the stochastic
this work that there is only one such point. We will also
SVIR model, the phenomenon of resonant response to
assume that Eq. (3) has one fixed point x in the hy-
vaccination. Section V contains concluding remarks. S
perplane x = 0. The state x is stable with respect to
E S
all variables except x . We call it the disease extinction
E
state. If x >0 (there is a nonzero number of infected),
II. THE DISEASE EXTINCTION RATE E
the deterministic trajectoryleavesthe vicinity ofx and
S
approaches the endemic state x .
A
A. The model of population dynamics Due to demographic noise the endemic state is ac-
tually metastable. A rare large fluctuation ultimately
We consider stochastic disease dynamics in a well- drivesthe populationinto a disease-freestate. The most
mixed population which includes infected (I) and sus- probablefluctuationbringsthe systemto the fixedpoint
ceptible (S) individuals and possibly other population x [13, 14]. The probability of such a fluctuation per
S
groups such as recovered or vaccinated. The system unit time, i.e., the disease extinction rate W , is given
e
state is described by a vector X = (S,I,...) with in- by the probability current to x , similarly to the prob-
S
teger components equal to the sizes of different popula- lem of escape from a metastable state [22]. For time-
tion groups. Along with X it is convenient to consider a independent W(X,r) this current is quasistationary for
3
timest ≪t≪W−1,wheret isthecharacteristicrelax- extinction rate, W ∝ exp(−Q), at N ≫ 1. This bar-
r e r e
ationtimeforthenoise-freemotiondescribedbyEq.(3). rier is entropic in nature, as it results from an unlikely
We note that, even though the state x is a saddle sequence of elementary transitions. It can be found by
S
point in the mean-field approximation, it differs from either solving the mean first passage time problem for
the saddle-point states encountered in the problem of reachingx (t)[7–9]orbycalculatingthetailofthequasi-
S
switchingbetweenmetastablestatesofreactionsystems. stationary probability distribution P(X,t) for x close to
In the case of interstate switching, the rates of elemen- x (t)[11,13,14]. Herewechoosethelatterstrategyand
S
tary transitions W(X,r) in the unstable direction are determine, to the leading order in N, the logarithm of
nonzero, and ultimately fluctuations drive the system the distribution tail. We seek the solution of Eq. (1) in
away from the saddle point. In the case of extinction, eikonal form, P(X,t) = exp[−Ns(x,t)] [28–30]. In the
fluctuations around x occur only in the extinction hy- limit of large N, from Eq. (1) we obtain the following
S
perplane, whereas the probability of exiting this hyper- equation for s(x,t):
plane is zero.
∂ s=−H(x,∂ s,t), (4)
If the system has, in the mean-field approximation, t x
morethanonesteadystate awayfromthe extinctionhy- H(x,p,t)= w(x,r,t)[exp(pr)−1].
perplane, extinction can go in steps: from the endemic Xr
state to another steady state and, ultimately to the ex- Here,wehavetakenintoaccountthat,typically,|r|≪N,
tinction hyperplane. In particular, if the only additional andW(X,r,t)dependsonXpolynomially,whereasP is
steady state is a saddle point at the boundary between exponential in X. Therefore we expanded P(X+r,t)≈
the basins of attraction of xA and xS, the problem of P(X,t)exp(−r∂xS)andreplaced,totheleadingorderin
extinction can be reduced to the problem of escape over 1/N, w(x−N−1r,r,t) by w(x,r,t).
this saddle point [23]. Equation (4) has the form of the Hamilton-Jacobi
equation for an auxiliary Hamiltonian system with
Hamiltonian H(x,p,t); s(x,t) is the action of this sys-
2. Periodically modulated systems tem. The Hamilton equations of motion are
x˙ = rw(x,r,t)epr,
The above picture can be extended to the case where r
X
the transition rates are periodic functions of time, p˙ =− ∂ w(x,r,t)(epr−1). (5)
x
W(X,r,t + T) = W(X,r,t). Periodicity of some of r
X
the rates in time is a natural way of modeling sea- Thesetrajectoriesdetermine,inturn,themostprobable,
sonal variations of epidemics [18]. The attracting so- or optimal,trajectoriesthat the system followsin a fluc-
lution of Eqs. (3), which describes the endemic state, tuation to a given state x at time t. We will calculate
is no longer stationary. We will assume that this solu- action s(x,t) using these trajectories and thus find the
tion, xA(t), is periodic in time with the same period T, exponent in the distribution P(X,t).
x (t+T) = x (t). The asymptotic disease extinction
A A
state x (t) is also periodic in time; it lies in the hyper-
S
plane x =0. 2. Boundary conditions for the optimal extinction
E
Themostprobablefluctuationwhichcausesextinction trajectory
of the disease corresponds to a transition from x (t) to
A
xS(t) [20]. An important characteristic of such a tran- To find the boundary conditions for Hamiltonian tra-
sition is the period-averaged disease extinction rate We. jectories (5), we note that the quasi-stationary distribu-
It can be introducedif the modulationperiod T ≪W−1 tion P(X,t) has a Gaussian maximum at X (t). This
e A
and,in addition, tr ≪We−1. In this case,fortime t such means that, close to attractor xA(t), action s(x,t) is
that tr,T ≪ t ≪ We−1, a quasi-stationary time-periodic quadratic in |x− xA)| for stationary systems, whereas
probabilitydistributionisformed,centeredatxA(t). The for periodically modulated systems s(x,t) = s(x,t+T)
probabilitycurrentfromxA(t)toxS(t)isalsoperiodicin is quadratic in the distance from trajectory xA(t) [31].
time, and the period-averagedvalue of this current gives On the Hamiltonian trajectories that give such action,
We [20], in a direct analogy with the problem of switch- the momentum p ≡ ∂xs → 0 for x → xA(t), and since
ing between metastable states in noise-drivendynamical x = x (t),p = 0 is a fixed point (a periodic trajectory)
A
systems [24–27]. of the Hamiltonian dynamics, the trajectoriesof interest
start at t→−∞,
t
B. Eikonal approximation
s(x,t)= dtL(x˙,x,t), (6)
Z
−∞
1. Equations of motion L(x˙,x,t)= w(x,r,t)[(pr−1)epr+1].
r
X
We will be interested in evaluating the disease extinc- In the Lagrangian L, Eq. (6), p should be expressed in
tion barrier Q which gives the exponent in the disease terms of x˙,x using Eq. (5). Since w ≥0, we have L≥0.
4
Extinction barrier Q is determined by Ns(x,t) for x 3. The t→∞ value of the momentum on the Hamilton
intheextinctionhyperplane,x =0. Inthe spiritofthe trajectory
E
eikonalapproximation,wehaveto findsuch(x,t) inthis
hyperplane that s(x,t) be minimal. The minimum de-
The momentum component p is generically nonzero,
terminestheboundaryconditionsfortheoptimalHamil- E
as found for stationary systems [9, 13, 14, 32]. For peri-
tonian trajectory of extinction, (x (t),p (t)). The
opt opt odically modulated systems, one can show that p 6= 0
condition that s(x,t) is minimal with respect to x E
i6=E by extending the arguments presented in Ref. [13, 32].
on the extinction hyperplane means that p = ∂ s → 0
i xi This amounts to showing that the stable manifold of the
for i 6= E as the trajectory (x (t),p (t)) approaches
opt opt periodic orbit(x (t),p=0)lies entirely inthe invariant
thehyperplane. Theminimumofs(x,t)withrespecttot S
hyperplane x = 0, p = 0, and, as a consequence,
withintheperiodofmodulationisreachedifH(x,p,t)→ E i6=E
does not intersect the unstable manifold of the periodic
0 as the trajectory approaches the hyperplane.
orbit (x (t),p = 0) . Such intersection is necessary in
A consequence of conditions H(x,p,t) → 0 and A
ordertohaveaheteroclinictrajectorythatwouldgofrom
p →0 is that the momentum component p remains
i6=E E (x (t),p=0) to (x (t),p=0).
bounded on trajectory (x (t),p (t)). Indeed, near A S
opt opt
the extinction hyperplane, xE ≪1, we have The hyperplane xE =0, pi6=E =0 is formedby trajec-
tories
x˙ = r w(x,r,t)epr
E E
r
X
≈xE rrE[∂w(x,r,t)/∂xE]xE=0epErE. (7) x˙i6=E = r[w(x,r,t)]xE=0ri, (9)
X X
Here,weassumedthatw(x,r,t)isnonsingularatx →0 p˙ = − [∂ w(x,r,t)] (epErE −1).
E E r xE xE=0
and, since w = 0 for x = 0 and r 6= 0 [cf. Eq. (2)], X
E E
we expandedw inx to the lowestorder. Letus assume
E
nowthat|p |→∞forx →0. Thenonlythetermwith The invariance of this hyperplane is a consequence of
E E
maximal −rE ≡ −rEm should be kept in the sum over Eq.(2),whichleadstop˙i6=E =0andx˙E =0forpi6=E =0
rE in Eq.(7); it is also clearthat pE shouldbe negative, and xE =0.
otherwise the trajectory would not approach xE = 0. Toprovethatthestablemanifoldof(xS(t),p=0)lies
FromtheHamiltonequationforpE andEq.(7)itfollows entirelyintheinvarianthyperplanexE =0,pi6=E =0,we
that dpE/dxE ≈ −1/xErEm. This relation, along with first show that the trajectories, which are described by
the explicit form of the Hamiltonian H, show that, if pE Eq.(9) andwhich startclose to the state (xS(t),p=0),
were diverging for xE → 0, the Hamiltonian would not approach this state for t → ∞. Then, since the dimen-
become equal to zero but would remain ≈∂w/∂xE with sion of the hyperplane xE = 0, pi6=E = 0 is equal to
thederivativecalculatedforxE =0andrE =rEm. This the dimension of the stable manifold of (xS(t),p = 0),
contradiction shows that the assumption |pE| → ∞ is we conclude that the stable manifold indeed lies in the
wrong, pE remains limited for xE →0. hyperplane.
Equation(7) shows that x → 0 exponentially as
E Equations (9) for x are the mean-field equations
t → ∞. As x approaches zero, variables x are i6=E
E i6=E in the extinction hyperplane x = 0, cf. Eqs. (3), and
approaching the equilibrium position in the hyperplane E
therefore x → (x (t)) for t → ∞. Linearization of
x = 0. This happens because p → 0 and the dy- i S i
E i6=E Eq. (9) for p about (x (t),p=0) gives
namics of x in the hyperplane is described by the E S
i6=E
mean-field equations, Eq. (3). Therefore,
∞ p˙ = − [∂ w(x,r,t)] p r . (10)
Q=Nsext, sext = dtL(x˙,x,t), (8) E Xr xE xS(t) E E
Z
−∞
x(t)→xS(t), p(t)→pS(t)fort→∞. We compare this equation with the mean-field equa-
tion for x near x (t). The latter has the form x˙ =
Function p (t) is periodic in time, with p = 0 and E S E
S i6=E
x r [∂ w(x,r,t)] . Since the state x (t) is
with hitherto unknown component pE(t), which is dis- E r E xE xS(t) S
cussed below. unsPtable in xE direction in the mean-field approxima-
tion, from Eq. (10) p˙ /p < 0. Therefore, all trajecto-
The optimal trajectory (x (t),p (t)) is a hetero- E E
opt opt
ries on the hyperplane x = 0, p = 0 close to the
clinic Hamiltoniantrajectorythatgoesfromperiodic or- E i6=E
state (x (t),p = 0) approach this state asymptotically
bit (x (t),p = 0) to periodic orbit (x (t),p ), and the S
A S S
ast→∞,andthusthe stablemanifoldof(x (t),p=0)
action for extinction s is calculated along this trajec- S
ext
lies in this hyperplane.
tory. Thetrajectoryx (t)istheoptimalpathtodisease
opt
extinction: itdescribesthemostprobablesequenceofel- From the above analysis one concludes that there
ementarytransitionsleadingtoextinction. Wenotethat, are no Hamiltonian trajectories that would go from
in periodically modulated systems, there is one optimal (x (t),p = 0) to x (t),p = 0 . Therefore the opti-
A S
path per period, whereas in stationary systems trajecto- mal trajectory lead(cid:0)ing to extincti(cid:1)on should go to a state
ries (x (t),p (t)) are time-translation invariant. x (t),[p (t)] with [p (t)] 6=0.
opt opt S S E S E
(cid:0) (cid:1)
5
III. OPTIMAL VACCINATION as
T
A. Constraint on the vaccination protocol T−1 dtξ0(t) W(1)(X,rξ,t) P(X,r,t)
Z0 XX (cid:12)(cid:12) (cid:12)(cid:12)
=NΞ. (cid:12) (cid:12) (11)
Vaccination increases the number of individuals who
are at least temporarily immune to the disease. It thus
Here, Ξ is the average vaccination rate rescaled by the
reducesthepoolofsusceptibleindividualsandultimately
characteristic population size N. The constraint is well-
leadstoareductionofthenumberofinfected. Whenthe defined for t ≪ t ≪ W−1, where P(X,r,t + T) ≈
r e
available amount of vaccine is small, so that the disease
P(X,r,t). Since for N ≫ 1 the population distribution
extinction still requires a large fluctuation, the goal of
sharply peaks at the endemic state X (t), the sum over
A
vaccinationis to reduce the disease extinction barrierQ. X in Eq. (11) can be replacedby W(1)(X (t),r ,t) , in
A ξ
Anoutcomeofvaccinationisoftenmodeledasthecre- the leading order in 1/N. (cid:12) (cid:12)
(cid:12) (cid:12)
ation of a sub-population of vaccinated individuals out In the presence of vaccination, one can still seek a so-
of susceptibles. The corresponding elementary transi- lution of the master equation in the eikonal form. The
tion rate is W(X,r) = ξ (t)X for r = −1, r = 1 exponent Q in the extinction rate is again given by the
0 S S V
and r = 0, where subscripts V and S refer to action of an auxiliary Hamiltonian system, Eq. (8). The
i6=S,V
vaccinated and susceptible individuals, respectively, and Hamiltonian now has the form
ξ (t) is the control field that characterizes the vaccina-
0
tion (subscript S should not be confused with subscript H(x,p)=H(0)(x,p)+ξ (t)H(1)(x,p),
0
S used to indicate the extinction state). Another model H(0)(x,p)= w(0)(x,r,t)(epr−1), (12)
is vaccination of newly arrived susceptibles, which leads r
X
to an effective reduction of the arrival rate µN. In this H(1)(x,p)=w(1)(x,r ,t)(eprξ −1).
ξ
model, the elementary transition rate for the arrival is
W(X,r) = N[µ−ξ0(t)] for rS = 1 and ri6=S = 0, with Ourgoalistofindthe optimalformofξ0(t)whichwould
ξ0(t)N being the change in the arrival rate due to vacci- minimizethediseaseextinctionbarrierQsubjecttocon-
nation. straint (11). Since w(1)(x (t),r ,t) is a known periodic
A ξ
function of time, we can equivalently search for the op-
We will consider a general model where vaccination
timal vaccination rate ξ(t) ≡ ξ (t) w(1)(x (t),r ,t) . It
modifies the rate ofanelementarytransitionofa certain 0 A ξ
type; the change of the population in the corresponding minimizes the functional (cid:12) (cid:12)
(cid:12) (cid:12)
transition is r . The field ξ (t) characterizing the vacci-
ξ 0 T
nation is assumed to be weak. The affected rate has the s˜ [ξ(t)]=s [ξ(t)]+λT−1 [ξ(t)−Ξ]dt,(13)
ext ext
form W(X,rξ,t) = W(0)(X,rξ,t)+ξ0(t)W(1)(X,rξ,t), Z0
with W(0) being the rate without vaccination. The vac- ξ(t)=ξ(t+T)=ξ (t) w(1)(x (t),r ,t) ≥0,
0 A ξ
cinationeitherincreasesordecreasestherate,asfortran- (cid:12) (cid:12)
(cid:12) (cid:12)
sitions from susceptibles to vaccinatedor for vaccination (cid:12) (cid:12)
whereλistheLagrangemultiplier. Thefunctionals [ξ]
ofnewlyarrivedsusceptibles,respectively. Therefore,we ext
is given by Eq. (8), where the Lagrangian corresponds
will assume without loss of generality that ξ (t)≥0 and
0
to the Hamiltonian (12) and depends on the vaccination
thatW(1)(X,r ,t)iseitherpositiveornegative. Wecon-
ξ
rate ξ(t).
sidermodelsinwhichthenumberofsusceptibleschanges
by 1 in an elementary transition associated with vacci-
nation, (r ) = ±1. We note that the analysis can be
ξ S
immediately extended to describe other processes, like B. Vaccination protocol for stationary systems
modification of the infection rates [3] or recovery accel-
eration by administrating medicine.
1. Double optimization problem
Itshouldbenotedthatthevaccinationmodeladopted
in this work is probabilistic by nature. An alternative Fora low averagevaccinationrate Ξit suffices to keep
is wherevaccinationis done ina pre-determinedfashion, intheactions [ξ(t)]onlytheleading-orderterminξ(t).
ext
when a certain number of individuals are vaccinated per Since Hamiltonian (12) is linear in ξ, this term is linear
unittimeatagiventime. Theanalysisofsuchdetermin- in ξ, too. In the spirit of the standard perturbation the-
istic vaccination lies beyond the scope of this paper. ory for Hamiltonian systems [33], the change in the ac-
tion,causedbythesmallperturbation,canbecalculated
We will assume that vaccination is periodic, ξ (t) =
0 alongtheunperturbedtrajectory(x(0)(t),p(0)(t))ofthe
ξ0(t+T), and that the amount of vaccine available per opt opt
period T is limited. We model this limitation as a con- Hamiltonian H(0), which describes the optimal path of
straint on the ensemble-averaged number of individuals disease extinction in the absence of vaccination.
vaccinated per period T. The constraint can be written Westartwiththecaseofsystems,whicharestationary
6
in the absence of vaccination. For such systems simple explicit form
sext[ξ(t)]=s(e0x)t+s(e1x)t[ξ(t)], (14) sext =mins˜ext =s(e0x)t+s(e1x)t,
s(1)[ξ(t)]=min ∞ dtχ(t−t )ξ(t), s(e1x)t =ΞT min χT(t). (17)
ext t0 Z−∞ 0 0≤t<T
−1 Alternatively, this expression can be rewritten in terms
χ(t)=−H(1) x(0)(t),p(0)(t) w(1)(x ,r ) .
opt opt A ξ oftheFouriertransformofthelogarithmicsusceptibility:
(cid:12) (cid:12)
(cid:0) (cid:1)(cid:12) (cid:12)
(cid:12) (cid:12)
Thequantityχ(t)iscalledlogarithmicsusceptibility[13, s(1) =Ξ min χ˜(nΩ)exp[inΩt], (18)
ext
20, 34, 35]; it gives the change of the logarithm of the t n
X
∞
extinctionrate,whichislinearinthevaccinationratefor
χ˜(ω)= dtχ(t)exp(iωt),
moderately low vaccination rate. Z
−∞
The minimization over t in Eq. (14) accounts for lift-
0
ing the time-translational invariance of the optimal ex- where Ω=2π/T is the cyclic frequency of vaccination.
tinction paths mentioned earlier. For ξ(t) ≡ 0, extinc- Weareinterestedinthesolutionforwhichs(1) isnega-
ext
tion can occur at any time (tr ≪ t ≪ W−1) with rate tive,whichrequiresminχT(t)<0. Onlyinthis casewill
W . Periodicvaccinationsynchronizesextinctionevents; vaccination reduce the barrier for disease extinction and
e
it periodically modulates the extinction rate, and the increase the disease extinction rate. The barrier reduc-
modulation is exponentially strong for N|s(1)| ≫ 1 (see tion due to vaccination, Q(1) = Ns(1) ∝ NΞ, becomes
ext ext
below). Formally, in a modulated system there is only large for N ≫1 even if the averagevaccinationrate Ξ is
one optimal extinction path per period, as explained in small. Thus, for not too small vaccination rates, where
Sec. II, which is here reflected in the minimization over the eikonal approximation is valid [20, 34], the effect of
t . This optimal path minimizes the disease extinction vaccinationonthediseaseextinctionrateisexponentially
0
barrier Q = Ns [19, 20, 34, 35]. Equation (14) is strong.
ext
closely related to the Mel’nikov theorem for dynamical Theexpressionfortheactionchanges(1),Eq.(17),can
ext
systems [20, 36]. be also obtained in a more intuitive way. Indeed, since
The constraint for minimizing the action over ξ(t) in ξ(t) is non-negative, it follows from Eq. (15) that
Eq. (13) has a form of an integral over the vaccination
period T. It is therefore convenient to write action s(1) (1) T
also in the form of such an integral, ext sext[ξ(t)]≥mtinχT(t)Z dtξ(t)=ΞT mtinχT(t). (19)
0
s(e1x)t[ξ(t)]=mt0inZ0T dtξ(t)χT(t−t0), TnisThme).imnFiimonriamml.auIlmlny,fiastcmptir,noivtiisidstehtdheebiynospξtt(aitmn)ca=eloΞpfaTtthiPmtehnawδt(htae−drjeutmχstiTsn(t+to)
∞
thevaccinationpulsessoastoincreasetheprobabilityof
χ (t)= χ(t+nT). (15)
T disease extinction. This provides the mechanism of syn-
n=X−∞
chronization by vaccination. Equation (19) immediately
The function χT(t) is obtained by superimposing the leads to Eqs. (16) and (17) with tλ replaced by tmin.
parts of χ(t) which differ by T. By construction, χ (t) In addition to the constraint on the average vacci-
T
is periodic in time t. nation rate, there may be an upper limit on the in-
stantaneous vaccination rate, which is imposed by con-
dition w(x,r) = w(0)(x,r ) + ξ(t)w(1)(x,r ) ≥ 0. In
ξ ξ
2. Temporal shape of optimal vaccination the case where w(1)(x,r ,t) < 0, as for vaccination of
ξ
newlyarrivedsusceptibles,thisconditionismetprovided
To find the optimal shape of vaccination rate ξ(t) we ξ0(t) ≤ ξ0m ≡ min{w(0)(x,rξ)/ w(1)(x,rξ) }. In this
first minimize the time integral in the variational prob- case the optimal vaccination proto(cid:12)col change(cid:12)s.
(cid:12) (cid:12)
lem Eqs. (13) – (15) with respect to ξ(t) for a given t . The new protocol can be found from the variational
0
Since ξ(t)≥0, it is convenient to perform the minimiza- problem(13)bychangingfromξ(t)≡ξ0(t) w(1)(xA,rξ)
tion with respect to ξ1/2(t) rather than ξ(t). The min- to an auxiliary function η(t) such that ξ0((cid:12)t) = ξ0m[1+(cid:12)
imization shows that ξ1/2(t) 6= 0 only for t = tλ, where η2(t)]−1, and then finding the minimum(cid:12)of s˜ext with(cid:12)
t is given by equation χ (t −t ) = −λ/T. From the respect to η(t). This choice satisfies the constraints
λ T λ 0
constraint on the period-averagedξ(t) we then have 0 ≤ ξ0(t) ≤ ξ0m. Variation with respect to η(t) shows
that s˜ has an extremum for η(t) = 0 or η(t) = ∞ for
ext
t6=t ,wheret isgivenbyequationχ (t −t )=−λ/T.
ξ(t)=ΞT δ(t−t +nT). (16) λ λ T λ 0
λ
The value of η(t) at the isolated instances t=t is arbi-
Xn λ
trary. Only regions where η(t) = 0, so that ξ (t) = ξ ,
0 0m
Substituting this expression into the functional s˜ and contribute to s˜ . Obviously, s˜ is minimal when
ext ext ext
minimizing with respect to t , we obtain the action in a ξ (t) = ξ for |t − t | ≤ ∆t/2, where t is the
0 0 0m min min
7
time whenχ (t)isminimaland∆tisdeterminedbythe 1. Limited lifetime of the vaccine
T
average vaccination rate Ξ. In other words, the vaccina-
tion rateξ(t) has the formof periodic rectangularpulses
The vaccination period T is naturally limited by the
of width ∆t, centered at t +nT, n = 0,±1,±2,....
min effectivelifetimeτv ofthevaccine. Thislifetimeisusually
The pulse width is
determined by the maximum storagetime of the vaccine
ΞT and/or by the mutation rates of the infectious agent. If
∆t= . (20)
ξ w(1)(x ,r ) τv is long, τv ≫ tr, vaccination can be made most ef-
0m A ξ
ficient by increasing the vaccination period up to ∼ τv.
(cid:12) (cid:12)
Since the vaccination (cid:12)rate is lim(cid:12)ited by the rate Indeed, as it follows from Eq.(17), s(1) ∝T in this case.
ext
of elementary transitions without vaccine, we have This result is easy to understand. Even though a de-
ξ0m w(1)(xA,r) . t−r1. Then for weak vaccination, crease of the vaccine pulse frequency Ω = 2π/T causes
ΞT (cid:12)≪ 1, from (cid:12)Eq. (20) ∆t ≪ tr. Therefore, χT(t) = a decrease of the prefactor in the disease extinction rate
(cid:12) (cid:12)
χT(tmin)during the pulse ofξ(t), tothe leadingorderin W , the exponential factor exp(−Ns(1)) in W increases
ΞT [we note thatχ (t)may varyona time scaleshorter e ext e
T sharply. Indeed, it can be seen from Eqs. (5), (12) and
thant ,seebelow;however,thistimescaleisalwayslong
r (14) that, as the system moves along the optimal path
compared to ∆t for sufficiently weak modulation]. The
to extinction, χ(t) is significant when the system is far
resulting change of the action is again given by Eq.(17).
fromthestationarystatesx andx . Thecharacteristic
A S
time scale of this motion is ∼ t . For T ≫ t we have
r r
min χ (t)≈min χ(t) and
C. Vaccination protocol for periodically modulated 0≤t≤T T t
systems
s(e1x)t =ΞT minχ(t), τv &T ≫tr. (21)
t
Optimalvaccinationinperiodicallymodulatedsystems
requiresaseparateconsideration. Here,thereisonlyone In the opposite limit of τv ≪ tr, and thus T ≪ tr, we
optimal extinction path per period T in the absence of have from Eq. (18)
vaccination. Whenthe averagevaccinationrateΞislow,
vtuarcbcinthaitsiopnatwhi.thTtohfierssatmoredpereriinodΞ,Tthweilllinoenalyriwneξaktleyrmpeirn- s(e1x)t =Ξχ˜(0), tr ≫τv &T, (22)
theactionstillhastheformofEq.(15),butwithoutmin-
In this case the vaccination-induced reduction of the ex-
imization overt . Since ξ(t)≥0, the minimum of action
0 tinction barrier is independent of the vaccination period
is still achieved for ξ(t) = ΞT δ(t−t +nT), but
n min and is determined by the zero-frequency component of
now tmin is uniquely determinePd by the strong modula- the logarithmic susceptibility.
tion. Inotherwords,themodulationuniquelydetermines
Aninterestingsituationmayoccurintheintermediate
the phase of the optimal vaccination pulses. The result-
ing expression for s(1) for the optimal vaccination pro- rangeτv ∼tr if,inthemean-fielddescription,thesystem
ext approaches the endemic state in an oscillatory manner.
tocol has the form of Eq. (17). If the vaccination pulses
In this case function χ(t) is also expected to oscillate.
are applied at a wrong time, i.e. if the phase difference
The oscillations are well-pronounced if their typical fre-
betweenthe vaccinationandthe modulationdiffers from quencyisω ≫t−1. ItisclearfromEq.(18)thatastrong
0 r
the optimal one, the vaccination will be not as efficient
effectondiseaseextinctioncanbeachievedbytuningthe
and may even be harmful: it may prolong the lifetime
vaccination frequency Ω = 2π/T or its overtones in res-
of the endemic state by increasing the diseaseextinction
onance with ω . An example of such a resonance will be
0
barrier Q.
discussed in Sec. V.
IV. DISEASE EXTINCTION BARRIER AS A
FUNCTION OF VACCINATION PERIOD 2. Limited vaccine accumulation
The vaccination-induced reduction of the disease ex- Adifferent situationoccursif the totalamountofvac-
tinction barrier Q(1) =Ns(1), as given by Eqs. (17), de- cine that can be accumulated is limited. This limitation
ext
pendsontheinterrelationbetweenthevaccinationperiod implies that ΞT ≤ M (note that M is the limit on the
T and the characteristic time scales of the logarithmic ensemble-averaged amount of the accumulated vaccine).
susceptibility χ(t). Function χ(t) may or may not oscil- Such a constraint is typical for live vaccines, as it may
late in time, but generally χ(t) is relatively large within bedangeroustostoretoomuchvaccineinthiscase. The
a time interval of the order of the relaxation time of the actual average vaccination rate in this case is now T-
system t [20, 34]. To revealsome qualitative features of dependent. We use the notation Ξ for this rate, with
r a
theeffectofvaccinationandinparticular,itsdependence Ξ =min(Ξ,M/T). This is Ξ that should be used now
a a
(1) (1)
on the vaccination period, we will consider s for two in Eqs. (21) and (22) for s in the limits T ≫ t and
ext ext r
types of constraint on this period. T ≪t , respectively.
r
8
The behavior of s(1) with varying vaccination period
ext
T depends on the form of the logarithmic susceptibil-
ity χ(t). Let us first consider the case where χ(t) has
(cid:0) (cid:9)
a single local minimum (at t = t ), and |χ(t)| mono-
∗
tonically decays to zero with increasing |t−t |. Here,
∗ (cid:5)(cid:15)(cid:6)(cid:7)(cid:8)
once the vaccine accumulation has reached saturation
(cid:11) (cid:12)(cid:13)(cid:9) (cid:14)
with increasing T (which happens for ΞT = M), func- (cid:4) (cid:10) (cid:1) (cid:2) (cid:3)
tion |s(1)|=M|minχ (t)| monotonically decreaseswith
ext T
further increase in T. Indeed, (cid:0) (cid:0) (cid:0) (cid:0)
d ∂
min χ (t)= χ(t −a T +nT)
T ∗ ∗ FIG. 1: The SVIR epidemic model with susceptible, vacci-
dT 0<t<T ∂T Xn nated, infected and recovered sub-populations. The arrows
dχ(t−a T +nT) indicateprocesses leadingtochangesof sub-populationsizes;
∗
= (n−a )>0,(23)
∗ thecorresponding rates are scaled perindividual.
(cid:18) dt (cid:19)
Xn t=t∗
where a gives the position of the minimum of χ (t)
∗ T For β > Γ ≡ γ + µ the SIR model possesses a sin-
over t and is given by equation ∂χ (t −a T)/∂a = 0;
T ∗ ∗ ∗ gle endemic state. This state corresponds, in the mean-
we choose 0 < a < 1 and take into account that
∗ field theory, to an attracting fixed point on the two-
min χ (t) < 0. In Eq. (23) we have used that,
0<t<T T dimensionalphase plane of susceptibles and infected. At
if χ(t) is minimal for t = t , then dχ/dt > 0 for t > t
∗ ∗ µ < 4(β − Γ)(Γ/β)2 this attracting point is a focus.
and dχ/dt < 0 for t < t . It follows from Eq. (23) that,
∗ The populations of susceptibles, infected and recovered
once the vaccine accumulation has reached saturation,
exhibit decaying oscillations in time as the system ap-
further increase in T will only reduce the effect of the
proaches the endemic state. It was found in Ref. [14]
vaccine. This result is understandable because, if T in-
that, in this parameter range, the populations oscillate
creasesbeyondM/Ξ,the actualaveragevaccinationrate
also on the optimal disease extinction path. These oscil-
Ξ decreases.
a lations are illustrated in Fig. 2.
A counterintuitive situation may occur if χ(t) is os-
We will now incorporate vaccination and introduce a
cillating. Here the inequality (23) may be violated. As
sub-population of vaccinated X = V. The vaccination
4
a result, the dependence of the effect of the vaccine on
is described by the transition rate W(X,r)=ξ (t)x for
0 1
T and, consequently, on the actual vaccination rate Ξ ,
a r = −1,r = 1,r = 0. The corresponding term in
1 4 i6=1,4
may be nonmonotonic. An example of this behavior is
the Hamiltonian Eq. (12) has the form ξ (t)H(1) with
discussed in the next section. 0
H(1)(x,p)=x ep4−p1 −1 . (24)
1
(cid:0) (cid:1)
V. RESONANCES IN THE STOCHASTIC SVIR Vaccinated individuals leave at the same rate µ as in-
MODEL
dividuals in other populations. For simplicity, we as-
sume that the immunity from the vaccination is never
Wenowapplysomeofourresultstoanimportantand lost. In this case fluctuations of the vaccinated pop-
widely usedstochasticepidemic model, the Susceptibles- ulation do not affect fluctuations of other populations,
Vaccinated-Infected-Recovered (SVIR) model. The and p ≡ 0 along the optimal extinction path. Then
4
modelissketchedinFig.1. Intheabsenceofvaccination, from Eq. (14), the logarithmic susceptibility is χ(t) =
ξ0(t)=0, the SVIR model reduces to the stochastic SIR x(0) (t)x−1 1−exp[−p(0) (t)] ,wherex(0) (t),p(0) (t)
modelwithpopulationturnover,whichwasoriginallyin- 1opt 1A(cid:16) 1opt (cid:17) 1opt 1opt
and x are calculated for the SIR model.
troduced to describe the spread of measles, mumps, and 1A
The Fourier spectrumofthe logarithmicsusceptibility
rubella, see [1, 3]. In the SIR model, susceptible in-
χ˜(ω)isplottedinFig.3(a). Itcorrespondstotheoptimal
dividuals are brought in, individuals in all population
extinction path shown in Fig. 2. As one can see, the
groups leave (for example, die), a susceptible individual
spectrum has a peak at the characteristic frequency of
can become infected upon contacting an infected indi-
oscillations of the system in the absence of vaccination
vidual, and an infected individual can recover. If we set
ω (for the chosen parameter values ω ≈5.2µ).
X = S,X = I, and X = R, the rates of the cor- 0 0
1 2 3
Wenowconsidertheeffectoftheresonantpeakinχ˜(ω)
responding processes are: (i) influx of the susceptibles,
on vaccination. The dependence of the scaled change of
W(X,r) = µN for r = 1,r = 0, (ii) leaving, with
the same rate for all1populati6=io1ns, W(X,r) = µXi for the disease extinction barrier s(e1x)t = Q(1)/N on vaccina-
r = −1,r = 0, (iii) infection, W(X,r) = βX X /N tion period T is shown in Fig. 3 (b). The solid line in
i j6=i 1 2
for r = −1,r = 1,r = 0, and (iv) recovery of the Fig. 3 (b) shows the behavior of s(1) where there is no
1 2 i6=1,2 ext
infected, W(X,r)= γX for r = −1,r =1,r = 0, limit on vaccine accumulation or, equivalently, for such
2 2 3 i6=2,3
Fig. (1). periodswherethelimitationdoesnotcomeintoplayand
9
theactualvaccinationrateΞa isindependentofT. Func- 2.5x 10−3 0
tion |s(e1x)t| ≡ −s(e1x)t is seen to be strongly nonmonotonic, χ~(ω′) (a) s′ext (b)
itdisplayspronouncedmaxima(whichcorrespondtothe 2
M′=1
minima of s(1)). They occur where the vaccination pe- −0.005
ext 1.5 ω′
riod T coincides with the multiples of the characteristic
0
period of the system motion without vaccination 2π/ω0. 1 2ω π′ M′=3
For limited vaccine accumulation M, the actual aver- −0.01 0
age vaccination rate depends on the vaccination period, 4 π
0.5 ω′
Ξa = min(Ξ,M/T). Beyond a certain value of T, the 0
ilneacdresatsoe oafcThainsgaecocofmthpeandieepdenbdyetnhceedoefcsr(e1a)seonΞaT..TRheis- 00 5 ω′ 10 −0.0150 2 4 6 T′ 8
ext
markably, |s(1)| still displays resonant peaks at 2πn/ω
ext 0 FIG. 3: (a) The Fouriertransform of thelogarithmic suscep-
with integer n. Their amplitude decreases with increas-
tibility intheSIRmodel. Theparameters are thesame asin
ingn. Theoccurrenceofthepeaksshowsthat,bytuning ′
Fig.2;therescaled frequencyisω =ω/µ. Thesusceptibility
the vaccination period, the effect of the vaccination can spectrumdisplaysasharppeakatthecharacteristicvibration
be resonantly enhanced; the resonance in this case is in frequencyω . (b)Thechangeofthescaled extinctionbarrier
0
the exponent of the disease extinction rate, and there- s′ = µs(1)/Ξ with vaccination period T. The solid line
ext ext
fore it is extremely strong. Counter-intuitively,since the shows s′ where there is no limit on vaccine accumulation,
ext
actualaveragevaccinationratedecreaseswithincreasing whereas thedashed linesrefer tothecase of limited accumu-
T, a strong enhancement of the vaccine can be achieved lated vaccine amount. The accumulation limit M is scaled
where this rate is decreased. For example, in Fig. 3 (b) bythesmall-T average vaccination rateΞ, M′ =µM/Ξ, and
′ ′
the maxima of |s(1)| for µM/Ξ = 1 and µM/Ξ = 3 are T =Tµ. Thelocationsoftheresonancesofsext areindepen-
ext dent of M.
achieved for T in the range where Ξ <Ξ.
a
0.08 fluctuation. We show that the optimal vaccination leads
x to anexponentially strongincrease ofthe diseaseextinc-
2
0.06 tion rate. This happens because the vaccine changes the
effective entropic barrier that needs to be overcome for
spontaneous extinction.
0.04
Wefindthattheoptimalvaccinationprotocolisaperi-
odicsequenceofδ-likepulses. Thisprotocolisessentially
0.02
model-independent, it only requires that the population
be spatially uniform. In stationary systems, the phase
0 of the pulses is irrelevant. In contrast, in periodically
0.4 0.5 0.6 0.7 0.8 0.9 1
x
1 modulated systems, like in the case of seasonally vary-
inginfection,itisnecessarytoappropriatelysynchronize
FIG.2: Themost probable trajectories in thestochastic SIR vaccinationpulseswiththemodulation. Moreover,ifthe
modelontheplaneofthescaled numbersofsusceptiblesand pulse phase is wrong, vaccination may hamper disease
infected, x1 = X1/N and x2 = X2/N, respectively. The extinction.
dashedlineshowsamean-fieldtrajectorytowardtheendemic
For fixed average vaccination rate, the effect of vacci-
state, and the solid line shows the most probable trajec-
nation in stationary systems increases with the increas-
tory followed during the fluctuation-induced disease extinc-
ing vaccination period. However, this increase is gener-
tion [14]. The plot refers to β/µ=80, and γ/µ=50.
ally nonmonotonic and the disease extinction rate can
display exponentially strong resonances. They occur if
the vaccination period coincides with the period of de-
caying oscillations of the population, which character-
VI. CONCLUSIONS ize the approach to the endemic state in the mean-field
(fluctuation-free) approximation. The resonances occur
We have developed a theory of optimal periodic vacci- also where the vaccination period coincides with a mul-
nation against an endemic disease for low average vacci- tiple of the dynamicalperiod. They are illustrated using
nation rate, as in the case where the vaccine is in short the well-known SVIR model of population dynamics.
supply, or short lived, or cannot be stored in the suf- Itturnsoutthat,counterintuitively,theeffectofvacci-
ficient amount. We assume that the vaccination rate is nationcanbe sometimes enhancedby reducingthe aver-
insufficientforeliminatingtheendemicstateandthusex- agevaccinationrate. This happens wherethe mean-field
terminating the disease by “brute force”. However, vac- dynamics is characterized by decaying oscillations and
cination can change the rate of disease extinction, which there is a constraint on the amount of vaccine that can
occurs spontaneously as a result of a comparatively rare be stored. In this case lowering the average vaccination
10
rate can allow one to tune the vaccination period in res- cinationratecannotbenegativeanditistheaveragevac-
onance with the system dynamics. cinationratethatisgiven. Theanalysiscanbeextended
The analysis is based on the master equation for the tootherproblemsofoptimalcontroloffluctuation-driven
population dynamics. We solve it in the eikonal approx- extinction with similar constraints.
imation and reduce the problems of the tail of the dis-
tribution andofthe extinctionto Hamiltoniandynamics
ofanauxiliarysystem. A generalformulationofthe cor-
Acknowledgments
responding Hamiltonian problem is obtained for period-
ically modulated systems. The optimal vaccination pro-
tocol is found using this formulation with account taken The workatMSU wassupported in partby the Army
of the constraint on the average vaccination rate. The ResearchOfficeandbyNSF GrantPHY-0555346. B.M.
feature ofthe problemthatmakesitdifferentfromother wassupportedbytheUS-IsraelBinationalScienceFoun-
problemsofoptimalcontrolofrareeventsisthatthevac- dation Grant 2008075.
[1] M. S. Bartlett, Stochastic Population Models in Ecology 011130 (2008).
and Epidemology (Wiley,New York,1960). [20] M. Assaf, A. Kamenev, and B. Meerson, Phys. Rev. E
[2] R. M. Anderson and R. May, Infectious diseases of hu- 78, 041123 (2008).
mans; dynamics and control (Oxford University Press, [21] D.RyvkineandM.I.Dykman,Phys.Rev.E73,061109
Oxford,1991). (2006).
[3] T.Andersson,H.&Britton,Stochastic EpidemicModels [22] H. Kramers, Physica (Utrecht)7, 284 (1940).
and Their Statistical Analysis, vol. 151 of Lecture Notes [23] M. Assaf and B. Meerson, Phys. Rev. E (in press);
in Statistics (Springer, New York,2000). arXiv:0907.0070.
[4] A. Scherer and A. McLean, British Med. Bull. 62, 187 [24] V.N.Smelyanskiy,M.I.Dykman,andB.Golding,Phys.
(2002). Rev. Lett.82, 3193 (1999).
[5] D. Robertson, B. Hahn, and P. Sharp, J. Mol. Evol. 40, [25] R. S. Maier and D. L. Stein, Phys. Rev. Lett. 86, 3942
249 (1995). (2001).
[6] A. Hay, V. Douglas, and Y. Lin, Philos. Trans. R. Soc. [26] J. Lehmann, P. Reimann, and P. Hanggi, Physical Re-
Lond.B Biol. Sci. 356, 1861 (2001). view E62, 6282 (2000).
[7] G.H.WeissandM.Dishon,Math.Biosci.11,261(1971). [27] M. I. Dykman and D. Ryvkine, Phys. Rev. Lett. 94,
[8] E. G. J. Leigh, J. Theor. Biology 90, 213 (1981). 070602 (2005).
[9] O. A. van Herwaarden and J. Grasman, J. Math. Biol. [28] R. Kubo,K. Matsuo, and K. Kitahara, J. Stat.Phys. 9,
33, 581 (1995). 51 (1973).
[10] I.N˚asell, J. Theor. Biol. 211, 11 (2001). [29] H. Gang, Phys.Rev.A 36, 5782 (1987).
[11] V. Elgart and A. Kamenev, Phys. Rev. E 70, 041106 [30] M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J.
(2004). Chem. Phys.100, 5735 (1994).
[12] C. R. Doering, K. V. Sargsyan, and L. M. Sander, Mul- [31] M. Dykman, X. L. Chu, and J. Ross, Phys. Rev. E 48,
tiscale Model. Simul. 3, 283 (2005). 1646 (1993).
[13] M. I. Dykman, I. B. Schwartz, and A. S. Landsman, [32] M. Khasin and M. I. Dykman, Phys. Rev. Lett. 103,
Phys.Rev.Lett. 101, 078101 (2008). 068101 (2009).
[14] A. Kamenev and B. Meerson, Phys. Rev. E 77, 061107 [33] L. D. Landau and E. M. Lifshitz, Mechanics (Elsevier,
(2008). Amsterdam, 2004), 3rd ed.
[15] B. Meerson and P.V. Sasorov, Phys. Rev. E 80, 041130 [34] V. N. Smelyanskiy,M. I. Dykman,H. Rabitz, and B. E.
(2009). Vugmeister, Phys. Rev.Lett. 79, 3113 (1997).
[16] V.N.SmelyanskiyandM. I.Dykman,Phys.Rev.E55, [35] V.N.Smelyanskiy,M.I.Dykman,H.Rabitz,B.E.Vug-
2516 (1997). meister, S. L. Bernasek, and A. B. Bocarsly, J. Chem.
[17] B. E. Vugmeister and H. Rabitz, Phys. Rev. E 55, 2522 Phys. 110, 11488 (1999).
(1997). [36] J. Guckenheimer and P. Holmes, Nonlinear Oscillators,
[18] N.C.Grassly andC.Fraser,Proc.R.Soc.Lond.B273, Dynamical Systems and Bifurcations of Vector Fields
2541 (2006). (Springer-Verlag, New York,1997).
[19] C. Escudero and J. A. Rodriguez, Phys. Rev. E 77,