A central challenge in disease ecology is developing effective management strategies to control harmful diseases and limit their damage to wildlife, human health, and ecosystem services. Using models of the dynamics of hosts and disease, this can be formulated as an economic optimization problem: how to maximize the benefits of a strategy, or how to minimize the costs of meeting a disease control goal.
Disease management may be considered an economic problem: Given limited resources, how to best manage disease to minimize impacts? Or, conversely, given management targets, how to best limit the cost of disease control? Research addressing such questions - economic epidemiology - is a growing field [@Brandeau2003; @REFS], though the literature to on the economics of wildlife disease is more limited. [@Horan2005; Bicknell1999; Fenichel2010]. The bioeconomics of invasive species management is more broadly explored and has many close parallels with the economics of disease management [@REF (Epanchin-Niell, Liebhold].
[TODO See additional citations in OPENQUALS for variety of methods applied to invasives]
Individual-based models (commonly called agent-based models) simulate of systems based on the behavior and properties of individuals, rather than the population as a whole. These can include network contact structure, spatial heterogeneity, individual variation, all of which affect disease transmission. Frequently such models can not be expressed as a closed-form set of equations at the population level without making strong assumptions about the distribution of individuals.
Importance of individual-based processes (network contact structure, host and pathogen variation) in disease dynamics. - IBMs work where population-level models do not capture appropriate dynamics or answer questions. [@Grimm2005]
Individual-based models have have been used in the management of disease in a "scenario" framework.
Agent-based models my include optimizing agents, rather than be optimized for management, in order to simulate the behavior of individuals.
One approach to optimization of IBMs has been to apply control techniques to reduced-form models. The optimization problem can be solved on these simplified ODEs. When the reduced-form model represents a mean-field approximation, it may perform poorly under when the underlying system is heterogenous [@Federico2013]. Another approach is to generate population-level models with non-mechanistic equations that approximate IBM behavior. Such methods can yield better results than mean-field models. However, the applicability of these adjusted models is limited the the range over which they were parameterized, and model parameters are difficult to interpret [Oremland 2015].
Optimal control, as practices by
@Pontryagin1967 described the necessary conditions to determine the optimal control path for a dynamical system. The application of these methods to dynamical ecological systems is detailed in @Clark1990. In brief, the optimal control path can be derived by solving maximization problem:
$$\max_{h >= 0} \int_{t=0}^T \pi(x, h, t) \, dt\, \text{ s.t. } \dot N, \dot P$$.
Various numerical methods can be used to find the optimal control path. Where $f(t)$ has a closed-form,
IBM dynamics can be difficult to capture as ODEs [@Brown2004]
Equation-free (EF), or multi-scale, modeling [@Armaou2004, but see @Kevrekidis2009a for review] is a method for capturing and analyzing population-scale dynamics of individual-based models while bypassing the deriving population-scale equations, which exist conceptually but often can not be derived. It has been used in a variety of applications, including physical chemistry, movement ecology and epidemiology [@REFS, Levin, @Williams2015]
In equation-free modeling, the population-scale system is simulated using repeated short-repeated bursts of the individual-based model. At each time step, a "lifting" operator maps the population-scale state to an ensemble of possible individual-scale states. Each of these instances of the individual-based model are simulated for a short burst of time, then are each reduced back to a population-level state with a "restricting" function. The results from all simulations are average to generate a mean result.
The system dynamics may be simulated by the lift-simulate-constrict cycle alone. However, it is computationally advantageous to use in conjunction with a projection step. In this approach, rather than simulating for an entire time-step, the simulation is run for a fraction of a time step to obtain an estimate the time derivative of the population-scale system. The population-scale state is then projected forward using this value in a differential equation solver using methods such forward Euler or Runge-Kutta. In addition these approximate derivatives may be used
For a stochastic system, the EF framework simulates the expected value of the simulation. Each micro-step consists of many simulation of
Originally developed for applications in physical chemistry, it has been applied to biological [@Erban2006], ecological [@Raghib2010a] and epidemiological systems [@Cisternas2004; @Williams2015].
The EF approach allows computationally efficient simulations of the expected value of population-level summaries of an individual based model.
The model is a host-pathogen system similar to the macroparasite model of @Anderson1978, where individuals can host multiple pathogens of the same type, and suffer increased mortality with greater infection load. The primary differences are the separation of birth and death processes and influx of infectious particles from outside the system.
In the mean-field version of the system, the host population $N$ increases via density-dependent reproduction ($rN(1-N/K)$) and decreases via a constant intrinsic mortality rate ($d$), and and additional mortality for each infection in a host ($\alpha$).
$$\frac{dN}{dt} &= r N (1 - N/K) - \alpha P - d N$$
The pathogen population $P$ grows with new infections that are density-dependent $\lambdaPN$, and decreases via pathogen mortality ($\mu P$), background host mortality ($dP$), and mortality due to the disease, which affects the most-infected hosts most $alpha(P + P/N)$. The form of of this term depends on the assumption of a random (Poisson) distribution of pathogens among hosts. Finally, additional parasites enter the system via external propagule pressure $\lambda_{ex}$. $\lambda_{ex}$ may be reduced by control effort $h(t)$, which reduces propugule pressure by a factor of $e^{h(t)}$.
$$\frac{dP}{dt} &= \lambda P N - \mu P - d P - \alpha P - \alpha P^2/N + e^{-h} N \lambda_{ex}$$
It assumption of (1) a large population size adequately represented by continuous variables, (2) a constant distribution of parasites among individuals, and (3) no stochastic processes.
The individual-based version of this represents the population as discrete counts of hosts and parasites. In the IBM, the state of the system is represented as a vector of individual hosts, each with a discrete number of infections. Births, deaths, new infections and loss of infections (recovery) in each individual $i$ with number of infections $j$, occur stochastically according to stochastic rates $r$:
$$\begin{aligned} r_{i, birth} &= r(1-N/K) \ r_{i, death} &= \alpha j + d\ r_{i, infection} &= \lambda * \sum_{i=1}^N j_i + lambda_{ex} e^{h(t)}, r_{i, recovery} &= \mu j \end{aligned}$$
This stochastic process occurs in continuous time, implemented via Gillespie's stochastic simulation algorithm (SSA) [@Gillespie1976].
Multi-scale simulation requires "lifting" and "restriction" functions, which map macro-scale system states to the micro-scale and vice-versa. In this case, I generate micro-level states from the macro state (lifting) by random sampling (with replacement) of $P$ individuals to distribute the individuals (a Poisson process). For the reverse (restriction), I simply sum the total host and parasite populations.
The management problem is to maximize the host population provides an ecosystem service of value $v$ per unit time, proportional to the host population size. The control $h$, which consists of reducing the arrival rate of new infections, has a cost $c$. Thus, the optimization problem is to maximize the net benefits over the course of a time period $T$, subject to the dynamics of the system.
$$\max_{h >= 0} \int_{t=0}^T \pi(x,h,t), dt\, \text{ s.t. } \frac{dx}{xt} = f(x, h, t) $$.
From this I derive the Hamiltonian equation, which must be optimized for all time points:
$$\begin{align} \mathcal{H} = &vN - ch + \ &\nu_1 \left( rN(1 - N/K) - \alpha P - d N \right) + \ &\nu_2 \left(\lambda P N - \mu P - d P - \alpha P - \alpha (P^2)/N + e^{-h} N \lambda_{ex} \right) \end{align}$$
Using the adjoint principal, I can derice the paths of the shadow values $\nu_1$ and $\nu_2$:
$$\begin{align} \frac{d\nu_1}{dt} &= -v - \nu_1 \left(r - d - 2r \frac{N}{K}\right) - \nu_2 \left(\lambda P + \alpha \frac{P^2}{N^2} + e^{-h} \lambda_{ex} \right) + \delta \nu_1 \ \frac{d\nu_2}{dt} &= \nu_1 \alpha - \nu_2 \left(\lambda N - \mu - d - \alpha - \frac{2 \alpha P}{N} \right) + \delta \nu_2 \end{align}$$
The IBM is an individual-based representation of the @Anderson1978 model of host-parisite dynamics.
ODE representation of the sytem makes assumptions of
ODE representation useful for equilibrium solutions, but not transients [@Adler1992; @Ross2015]
ODE system is reduced in the following ways:
This IBM is relatively simple. It lacks several features that are trivial to incorporate into an IBM but challenging to represent in closed-form solutions, such as individual variation in traits or heterogenity in the host contact structure.
Economic control problem:
[EQUATIONS]
For the equation-free simulation approach, the system dynamics are a "black box" which may be simulated but under which the derivatives $\dot N$ and $\dot P$ are unknown.
At each time step, I simulate the system using the lift-simulate-restrict cycle repeatedly in order optimize
The derivative of the shadow states can be derived via the adjoint principle:
$$\frac{d\nu}{dt} = - \frac{d\mathcal H}{dx} = -\frac{d\pi}{dx} - \nu \frac{df}{dx}$$
The profit function $\pi(\dot)$ is known, so $\frac{d\pi}{dx}$ can be determined analytically. To calculate $\frac{df}{dx}$,
While the equation-free framework can be used on "black box" simulators, taking advantage of some properties of the individual-based model can improve computational efficiency and performance of the method. In this case, I do so in two ways. First, as the IBM uses the Gillespie SSA, I use a single variable-time "Gillespie" micro-step rather than a fixed-period micro-step. This allows me to estimate $dx/dt$ without any micro-step offset. Second, I take advantage of the mapping between continuous population states at the macro-level and discrete counts at the micro-level. When lifting a continuous macro-level to create an ensemble of random micro-level states, some will have populations rounded down from the continuous level, and some rounded up. The difference in $x'$ between these populations is used to determine $x'/x$ without additional simulations.
For a single variable system $df / dx$, can be approximated if the second derivative of the state variable is calculated by the integrator. That is $$\frac{df(\dot)}{dx} / dx \approx \frac{\Delta f(\dot) / \Delta t}{\Delta x / \Delta t}$. However, for higher-dimensional systems, this approximation fails. Instead the full matrix of $$\frac{df(\dot)}{dx}$$ (the Jacobian matrix) is approximated by perturbing the simulation and measuring the difference in outcomes.
Simulate optimal path by solving for Hamiltonian maximum at each time step
Lifting and restriction functions translate between full distribution of infections and summary statistics of distribution
In the absence of control, the the disease invades the system rapidly. This suppresses the population via and increased overall mortality rate for hosts, and the system reaches a stable equilibrium.
In the ODE system, the internal equilibrium can be calculate as [TODO].
The individual-based, stochastic model results in figure. Shown are the mean path and the 95% confidence intervals for the system state through time,
EF version of IBM correctly produces results similar to overall IBM:
IBM and ODE have different responses to the the same controls.
EF optimization yields better results for the IBM than ODE-derived control
Complement to discrete approaches, ensure agreement with theory.
Still a model reduction approach, but reduction takes place at a different scale
Challenges
Stochasticity: The equation-free model produces an estimate of the IBM's mean behavior, and thus an optimal control path for the mean conditions. The primary utility of the method is not in stochastic control. [See a host of other methods for this: Talk to Carl/Jim]. Rather, its value lies in optimizing where the IBM behaviors are not easily expressed in closed form equations. Stochastic behavior makes the method computationally difficult. The variation in outcomes is amplified in the estimate of the derivatives and even more so in the estimation of of the Jacobian matrix.
Journals: Royal Society Interface? Ecological Economics? Methods in Ecol and Evo. AmNat for sepearate disease paper.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.