OPTIMAL CONTROL ON HYBRID ODE SYSTEMS WITH APPLICATION TO A TICK DISEASE MODEL

We are considering an optimal control problem for a type of hybrid system involving ordinary differential equations and a discrete time feature. One state variable has dynamics in only one season of the year and has a jump condition to obtain the initial condition for that corresponding season in the next year. The other state variable has continuous dynamics. Given a general objective functional, existence, necessary conditions and uniqueness for an optimal control are established. We apply our approach to a tick-transmitted disease model with age structure in which the tick dynamics changes seasonally while hosts have continuous dynamics. The goal is to maximize disease-free ticks and minimize infected ticks through an optimal control strategy of treatment with acaricide. Numerical examples are given to illustrate the results.

1. Introduction.In many economic, biological and physical scenarios, one needs to make reasonable choices at a number of decision making levels.While the specifics of the problem depend on the application context, the goal is to "optimize performance".The underlying system often has continuous dynamics, frequently governed by ordinary or partial differential equations.Decisions, such as resource allocation or control adjustment, may occur at discrete times.Thus these hybrid systems involve continuous and discrete components.The analysis needed for optimal control of such systems is a mixture of control techniques for infinite dimensional systems [10,18,32,34] and ideas from Pontryagin's Maximum Principle [27].
After developing the control analysis of a particular hybrid system, we will first be applying this approach to a model of a tick-transmitted disease-Lyme diseasewhich is a serious health problem affecting humans as well as domestic animals in many parts of the world.The infections are generally transmitted through a bite of an infected tick, and it appears that most infections are widely present in some wildlife species.Hence, an understanding of tick population dynamics and its interaction with hosts is essential to understand and control the disease [17].We use the model presented by Ghosh and Pugliese [13], which is a type of hybrid system with eleven ordinary differential equations.The dynamics of ticks (Ixodes ricinus in Trentino, Italy) changes seasonally and has some discrete components, while the hosts have continuous dynamics throughout the years.The goal is to maximize the disease-free ticks and minimize the infected ticks through an optimal control strategy of treatment with acaricide.In numerical examples, we also illustrate the case when only minimizing the infected tick and the cost of applying the acaricide.The most consistently effective method for reducing an abundance of ticks on residential properties is to spray or otherwise broadcast acaricides onto vegetation where the ticks live.Acaricides can be delivered directly to tick hosts to kill ticks on the animals.The use of a bait box to attract rodents and directly treat them with the acaricide fipronil kills ticks on the rodents.Another baited device to lure deer and treat them with acaricide has the potential to reduce tick populations over wide areas [15].
Ixodidae ticks have a two-year life cycle (Fig. 1).After hatching from eggs, they pass from one life stage to another by molting, after a blood meal.In temperate climates, the life cycle is strongly influenced by the seasonal rhythm [1,28].Eggs are laid by an adult female tick in the spring and hatch into larvae later in the summer.These larvae reach their peak activity in August.The larva then attaches itself to its host, begins feeding, and over a few days, engorges (swells up) with blood.Most larvae, after feeding, drop off their hosts and molt, or transforming into nymphs in the fall.The nymphs remain inactive throughout the winter and early spring.In May, nymphal activity begins.The nymph then latches on to its host and feed for four or five days, engorging with blood and swelling to many times its original size.Once engorged, the nymph drops off its host into the leaf litter and molts into an adult.These adults actively seek new hosts throughout the fall.Peak activity for adult ticks occurs in late October and early November.As winter closes in, adult ticks unsuccessful in finding hosts take cover under leaf litter or other surface vegetation, becoming inactive in temperatures below 40 degrees F. Adult female ticks that attach to deer, whether in the fall or spring, feed for approximately one week.Mating may take place on or off the host, and is required for the female's successful completion of the blood meal.The females then drop off the host, become gravid, lay their eggs underneath leaf litter in early spring, and die.Each female lays eggs which hatch later in the summer, beginning the two-year cycle anew [1].
Infection is transmitted from infected ticks to susceptible hosts, or vice versa from infected hosts to susceptible ticks, while a tick is feeding on a host.A larva feeding on an infected host may become, after molting, an infected nymph; analogously, a nymph feeding on an infected host may become an infected adult.In both cases, infection is assumed to last forever [13].On the other hand, we assume that a host, after a period of infection, will become immune and no longer capable of transmitting the infection.
There has been work on diseases in ticks from different perspectives.Sandberg et al. [31] used a matrix model to investigate the seasonally varying population densities of questing ticks.Awerbuch-Friedlander et al. [2] studied a nonlinear system of difference equations that models the three-stage life cycle of the deer tick over four seasons and showed that seasonality can increase the stability of the system.Mount et al. [24,25,26] formulated a series of computer simulations based on age-structured difference equations to examine the relationship between host density, tick density and the persistence of a tick-borne disease.Gaff and Gross [12] studied a tick-borne disease model incorporating nonconstant population sizes and spatial heterogeneity utilizing a system of differential equations that may be applied to a variety of spatial patches.By creating a set of patches that reflect an area of interest, Gaff [11] explored different control options to aid in choosing the most effective program prior to field application.Buskirk and Ostfeld [5] used a computer model to show that the density of ticks was more sensitive to the availability of hosts for juveniles than of hosts for adults; hence, modifying the average reservoir competence of hosts for juveniles, could be used to manipulate the risk of Lyme disease without actually changing the density of ticks.They also studied the importance of spatial heterogeneity by combining field sampling and modeling approaches [6].Busenberg and Cooke [4] studied vertically transmitted disease; i.e., the passing of the infection from parent to offspring.They presented a difference equation model for the transmission of rickettsia in ticks.They concentrate on one aspect of vertical transmission, which is the cumulative effect of some pathogens as they are vertically transmitted through several generations.Here we focus on the seasonal pattern of ticks in various life stages and their interactions with hosts and try to find a best strategy to control the disease.
In work related to optimal control of hybrid systems, Miller and Rubinovich [23] studied optimal impulse control problems with a restricted number of impulses.They used the method of discontinuous time change to derive the necessary and sufficient optimality conditions in the form of a maximum principle.Cassandras et al. [8] presented a hybrid system modeling framework (motivated from manufacturing environments) which combines the time-driven dynamics of various physical processes with the event-driven dynamics describing switches between the physical processes.We call attention to some papers in two special issues on "control of hybrid systems" in Automatica and IEEE Transactions of Automatic Control [16,21].Many of these papers deal with systems for which time is the underlying variable, entering in some states as continuous variables and in other states as discrete variables; there is some emphasis on stability and reachability issues [22,33].See the treatment of optimal control of systems with continuous time dynamics between jumps at discrete times in [7,19].
This paper is organized as follows: First we describe the tick-transmitted disease model.Then we set up the mathematical framework on a prototype model to deal with the optimal control of the hybrid ODE system, and derive the existence, necessary conditions, and uniqueness for the optimal control.Finally, we apply W. DING our approach to the tick-transmitted disease model and give numerical examples to illustrate the results.
2. Tick model with infection.As seen in the life cycle of ticks, there are two distinct periods in which ticks are active or dormant.So, we consider a simple model for tick dynamics with only two distinct time periods, called "summer" and "winter".To be precise, assume that larvae and nymphs that feed during one "summer" go through the molting stage but arrest their development in "winter" and emerge (as nymphs or adults, respectively) in the following "summer".On the other hand, assume that after the adult females feed and produce eggs, a proportion of the eggs hatch immediately, so that larvae emerge in the same "summer", while the rest arrest their development and larvae emerge in the following "summer".Finally, assume that larvae, nymphs, and adults die at the end of the "summer" in which they have emerged, if they have not succeeded in feeding.Of course, these assumptions are rather crude with respect to the complex interactions among climate factors, individual fat reserves, and feeding time; we believe, however, they capture some essential features of ticks' seasonal rhythm [13].
For the sake of simplicity, we disregard transmission between cofeeding ticks (i.e., direct transmission among ticks feeding close to each other).Also we assume that the infection does not affect either tick or host demography.Rodents have often been implicated as one of the important hosts for nymphs, while deer or other large mammals are preferred at the adult stage [3].Hence, for host 1 (rodents), we assume the total population is constant, and so is host 2 (deer).Otherwise, population might fluctuate for reasons other than interactions with ticks [13].Moreover, as in model [13], only host 1 (rodents) is assumed to become infected and divided into susceptible, infected, and recovered classes; host 2 (deer) is needed only to sustain the tick population which is not partitioned into those three classes, and host 2 remains as a constant.
Our state variables are as follows: where time t is the underlying variable.For host 1, we denote the total population as H 1 , which is taken as a constant; i.e., H 1 = H s n + H i n + H r n , but each class has its own dynamics.For host 2, we denote the total population as H 2 , which is also taken as a constant and is not divided into susceptible, infected, and recovered classes.
The assumptions mentioned above translate into the following model from Ghosh and Pugliese [13] supplemented with control coefficients.We first list the differential equations for ticks in the "summer": ( where • p is the probability of immediate development of tick larvae; • c is the average number of eggs produced per fed adult; • u n is the death rate of ticks of various stages in summer n, for n = 0, • • • , N , our bounded control variable.We control it by applying the acaricide to the ticks, the lower bound is the natural death rate of ticks and the upper bound depends on the effectiveness of the acaricide; • g z are the feeding rates in stage z, where z = L, N, A. As in Rosa et al. [30], the feeding rates g z are assumed to depend on host densities according to a saturation function, because of the extended feeding period; i.e., where , β z i are the rates of encounters between questing ticks in stage z, for z = L, N, A and hosts H i , i = 1, 2; σ z i are the detachment rates of ticks in stage z, z = L, N, A, feeding on host i, for i = 1, 2.
Note we assumed host densities to be constant; thus, feeding rates will also be constant.And the rate at which susceptible hosts become infected is , where q z is the probability of becoming infected for a host of type 1 by an infectious tick in stage z, with z = N, A.
This gives rise to the following equations for host 1: where a 1 is the birth rate for host 1 (rodents) and b 1 is the death rate for H 1 .We choose During "winter", ticks arrest their development, so we only have host 1 (rodents) dynamics, given by: (3)

W. DING
System (1) -( 3) needs to be complemented with initial conditions.Each class of host 1 has continuous dynamics; host 2 (deer) is constant, which is used in the feeding rate.We need to keep track of the fact that larvae and nymphs that have fed on infected hosts will emerge as infected nymphs and adults respectively.These initial conditions include some jumps, which are hybrid-type features.
• w is the probability of survival through winter for larvae that have delayed development; • m z are the molting rates in stage z with z = L, N, A, taken as constants; • T is the length of "summer".
We illustrate the initial condition for susceptible nymphs in the beginning of the (n + 1) th summer as follows: after molting, larvae will become susceptible nymphs, but since some larvae encounter an infected host of type 1, they will become infected nymphs.
Now we are ready to set up the optimal control problem.Given the control set We choose the lower bound to be the natural death rate of ticks and the upper bound to depend on the effectiveness of the acaricide.We want to maximize the following objective functional: where B, C, K are positive constants balancing the three parts of the objective functional.This means we want to maximize disease-free ticks and minimize infected ticks during N + 1 years while minimizing the cost of applying acaricide.We take a quadratic cost for simplicity; other formats could also be treated.If minimizing the infected ticks is the main goal, we can take C = 0 (done in numerical examples).
3. Framework of mathematical approach.Now we set up the framework to do optimal control of the hybrid ODE system.We need to develop the optimal control analysis for the seasonal pattern of tick life and the integral initial conditions because we cannot use Pontryagin's Maximum Principle [27] directly.
We consider an abstract model for the dynamics of two populations with two distinct seasons: "summer" and "winter".We consier N + 1 years and use subscript n to denote the population in year n, where n = 0, 1, )) has dynamics only in "summer" and has a jump condition to obtain the initial condition for next "summer"; population 2 (y(t) = (y 0 (t), y 1 (t), • • • , y N (t))) has continuous dynamics throughout the years.Then populations 1 and 2 have dynamics similar to ticks and hosts respectively.The corresponding state variables x(t), y(t) satisfy the state system: Initial conditions: where x 0 (0), y 0 (0) are given, and p 1 , p 2 , p 3 are constants.The initial condition for y means it has continuous dynamics, and 365 is the last day of the current year, which is also the first day of the following year.For the tick application, p 2 and p 3 type jump conditions are needed, but p 1 is not needed.Note p 1 is included here since it is a common condition in many other applications.The jump with p 1 means the information at the beginning of the next year depends on the information at the end of the previous year, a condition is often used in resource allocation in economics.In this section, our basic assumption is that are bounded and Lipschitz continuous for their arguments.So the solution to ( 7) -(9) exists and is a priori bounded and unique [20].
The goal is to maximize the objective functional: where 4. Existence of an optimal control.We use ideas from Fleming and Rishel [10] and Strauss [34] to prove the existence of an optimal control.THEOREM 4.1.In (7) − (10), ) is a concave function.Then there exists an optimal control u ∈ U M that maximizes J(u).
Proof.Since for any control u ∈ U M , m ≤ u n ≤ M, there exists C 1 such that the corresponding states satisfy Then sup u∈U M J(u) exists, and there exists a maximizing sequence Here we explicitly show the dependence of the states on the controls, and we note that time is still the underlying variable.Because of the uniform boundedness of , N , passing to the limit, we have Notice in the initial conditions (9), we can write So we obtain (x * n , y * n ) are solutions of ( 7) -( 8) associated with u * and satisfy the initial conditions (9).
Because k n (x n (t), y n (t), u n (t)) is concave, using upper semicontinuity with respect to weak convergence, Therefore u * is an optimal control that maximizes the objective functional J(u).

5.
Necessary conditions of optimal control.We need to differentiate the objective functional with respect to the control u.Denote Since x = x(u), y = y(u) are involved in J(u), we first state the appropriate differentiability of the mapping u −→ (x(u), y(u)).The proof of this "sensitivity" result follows from estimates on difference quotients x(u + v) − x(u) and y(u + v) − y(u) of the state system.
LEMMA 5.1.(Sensitivity) The mapping is differentiable in the following sense: there exist functions t ∈ [365n + T, 365(n + 1)], where with respect to y n .
To obtain the characterization of the optimal control, we need the following adjoint equations.We define λ THEOREM 5.2.Assume [H1], given an optimal control u * and the corresponding solutions x * , y * of the state system ( 7) -( 9), there exist adjoint variables λ 1 , λ 2 satisfying: ) and λ 2 solving: with Proof.Using ( 7) -( 8), we add three zero terms to the objective functional and then use integration by parts: We first consider boundary terms for population 1 in "summers". . . .
We reindex for n and obtain the boundary terms as x n (s)ds + p 3 365n+T 365n x n (s)y n (s)ds . ( Notice population 2 has continuous dynamics throughout the years, so the boundary terms for population 2 are simply: − λ 2n (365(n + 1))y n (365(n + 1)) + λ 2n (365n + T )y n (365n + T ) =λ 20 (0)y 0 (0) − λ 2N (365(N + 1))y N (365(N + 1)).(26) Then we rewrite the objective functional ( 21) using ( 25) and (26) as Let u * ∈ U M be an optimal control and (x * , y * ) be the corresponding optimal state solutions.Let be the corresponding solutions of the system ( 7) -( 8) for n = 0, • • • , N .We compute the directional derivative of the objective functional J(u * ) with respect to u * in the direction of h.Since J(u * ) is a maximum value, we have Using the chain rule, the sensitivity from Lemma 5.1; i.e., and and t ∈ [365N, 365N + T ], For each n = 0, • • • , N, on the set where u * n = m, we choose variation h n with support on this set and h n ≥ 0, which implies that This gives the desired characterization of the optimal control.6. Uniqueness of the optimal control.We show uniqueness of the optimal control by considering the convexity of the objective functional with respect to the control.THEOREM 6.1.Suppose in (7)− (8), for their arguments, their second partial derivatives are bounded, (k n ) uu < 0 with sufficiently large absolute value; then there is a unique optimal control.Proof.We have shown the existence of an optimal control; now we show uniqueness of the optimal control by showing for all u, v ∈ U M , 0 < < 1, g ( ) < 0, where g( . This implies the strict convexity of the following map: δ , as δ → 0 and the sensitivities ψ n , φ n satisfy: where the derivatives of f n , g 1n , g 2n are evaluated at ( From ( 34) and ( 35), we obtain where C 1 = C fn , C 2 = max{C fn , C g1,n }, and C fn , C g1,n are the bounds for the partial derivatives of f n (x n , y n , u n ) and g 1,n (x n , y n ) respectively, since we have Lipschtiz conditions.Gronwall's Inequality implies for 0 ≤ t ≤ T, For T ≤ t ≤ 365, we get a similar estimate from (36): where C 3 = (1+C g 2,n (365−T )e Cg 2,n (365−T ) ), and C g 2,n is a bound for the derivative of g 2,n (y n ).Since ψ 0 , φ 0 are bounded, we can use the bounds of To calculate g ( ), we need a second derivative of x, y with respect to the control u.As in Lemma 5.1, we have, where the derivatives of f n , g 1n , g 2n are evaluated at ( We can get the estimates for Then we are ready to calculate the derivatives of g.We have where the arguments for ( For the second derivative, we have (53) Note that for our objective functional which is quadratic in the control, we can explicitly solve the characterization for the optimal control: for n = 0, 7.2.Numerical results and discussion.Table 1 is the list of all parameters used in our numerical simulations.All the parameter values, which are considered to be reasonable for describing Ixodes ricinus tick populations in Trentino, Italy (see CEA Report 2000 for background information [29]), are from Ghosh and Pugliese [13], except w, C, B, m, M, m z , K. In all parameters time is measured in days, and host densities are per hectare.We take the probability of survival through winter for larvae that have delayed development w to be 0.5.For C, B, K, constants that balance the three parts in the objective functional, we take C = 1, B = 30 to indicate minimizing infected ticks is more important than maximizing the disease-free ticks.Other values can be investigated.Our goal is to see how the change in the cost of applying acaricide will affect the control strategy.We take the weight in the objective functional for the cost K to be 500, 5000, 50000.We choose m = 0.01 and M = 0.9 as the lower and upper bound for control u.The molting probability m z in [13] is taken as a given by 0.15e −0.008x , here we take it as a constant (m = 0.1) to derive the adjoint system using Theorem 5.2.
The optimality system consists of the state ODEs and the adjoint ODEs coupled with the optimal control characterization.We solve the optimality system of hybrid ODEs by the following forward-backward sweep method [14].Our initial conditions are (L 0 , N s 0 , N i 0 , A s 0 , A i 0 ) = (150, 200, 175, 20, 10) and (H s 0 , H i 0 , H r 0 ) = (20,5,5), which are the steady-state solutions from [13].We did the simulation for two years (730 days), because ixodidae ticks have a two-year life cycle, and the time step within each Rugga-Kutta iteration is 0.01.To see the figures clearly, our x-axis is [-10, 730] instead of [0, 730].
One might have two different goals: maximizing disease-free ticks and minimizing infected ticks (we took C = 1, B = 30 in this case); or only minimizing the infected ticks (we took C = 0, B = 30 in this case).Figure 2 is the optimal control for C = 1 and C = 0 respectively.We find in the C = 0 case we need to apply more acaricide (u = 0.61) in the beginning of the second summer with K = 500, while when C = 1, we only need u = 0.39.
Also we see the seasonal pattern of applying the acaricide, which is indicated by the seasonal pattern of ticks.We see the dramatically changing strategy of applying the acaricide.We apply it heavily for a short period of time in the beginning of the first summer, then apply the minimum amount (lower bound of the optimal control) until the end of summer; in the beginning of the second summer, we apply much less acaricide, then quickly reduce to the lower bound.We vary K to see how the change in the cost of acaricide affects the control u.We can see that increasing the cost will decrease the optimal control u, because it is more expensive to apply acaricide at higher cost.
We compare numerical results without and with control in Figures 3 -10.Without any control, we can see the tick's population in every life stage increases rapidly.We can see the optimal control strategy can significantly reduce the infected ticks density in every life stage and the density of infected host 1, and it also can increase the density of susceptible host 1.We observe a sharp change in the slopes of the tick density curves when the control hits its lower bound as a result of the control strategy.
Figure 3 is the graph for density of larvae.We can see when the cost is relatively high, the larvae density is also relatively high, because it is expensive to kill the ticks.The "hump" in the first summer is due to the adults producing new eggs emerging to larvae.
For Figures 4 and 5 we use log scale for y-axis because the differences in densities of the nymphs between the first and second summer are large.We can see that decreasing the cost will decrease susceptible nymphs, because it is cheaper to apply control.We can see the nymphs are killed quickly with the optimal control strategy.And we observe the drastic change in the slopes of the density curves within each summer when the control hits its lower bound.Figure 4 is for susceptible nymphs.For Figure 5, similar situations occur.As the cost of applying acaricide increases, we can see higher density of infected nymphs, because it is harder to remove them.Also we see the drastic change in the densities of infected nymphs during each summer.The infected ticks survive for only a short period.
Figures 6 and 7 also use a logarithmic scale for the y-axis.Figure 6 is for susceptible adult ticks.We can see that decreasing the cost will decrease both susceptible and infected adult ticks because it is easier to apply control.Figure 7 shows infected adult ticks.We see the changing behavior of the density of adult ticks.The graphs of infected host 1 are in Figure 9.We can see decreasing the cost will decrease the infected host-1 density, because we have fewer nymphs and adult ticks (see H i eqn).Figure 10 shows the recovered host 1.We can see decreasing cost  We developed the necessary conditions for an optimal control and corresponding states and adjoints for an unusual ODE system with a hybrid feature.The specific hybrid structure was motivated by a tick model from [13].We successfully applied our control results to determine optimal acaricide levels in this tick model and illustrated our results with numerical simulations.The optimal controls and corresponding states responded as expected to changes in the cost of acaricide, meaning higher cost coefficient gave lower level of control.We found that applying the acaricide treatment more heavily at the beginning of the first summer for a short period of time, then quickly reducing to the lower bound, applying much less in the beginning of the second summer then reducing to the lower bound, was optimal.This indicates the density of ticks in every stage can be decreased dramatically after a short period.And we see the drastic change in the densities of ticks in every life stage in each summer resulting from the optimal control strategy.
Notice we have only host 1 (rodents) dynamics in the model.This is the first step in understanding the interactions between the host and ticks to find the optimal strategy to control Lyme disease.Later we can include deer dynamics as well.
Our goal is to maximize disease-free ticks and minimize infected ticks considering cost to apply the acaricide.We can thus have different values for the weight coefficients for each component.In particular, we include graphs for only minimizing infected ticks (take C = 0) which required more control efforts in the beginning of the second summer when K = 500 and can reduce tick population more than in the C = 1 case.
The term "hybrid" has wide applications these days and can be used for a variety of systems that have components with distinct features.The existence, uniqueness, and robustness properties depend on the type of features involved.Our particular "hybrid system" has continuous dynamics during certain seasons and then has jump dynamics to start the next corresponding season.
In our numerical simulation, the optimality system is robust in the sense that it is not overly sensitive to the variation of the parameter values we have tested.If the functions on the right-hand side of the Hybrid ODE system are bounded and Lipschitz continuous for their arguments, then the solution exists and is a priori bounded and unique [20].

Table 1 .
Parameter Values

1 .
Initialization step: Choose initial guesses for ticks, hosts and control u.Update the control by entering new ticks, hosts values and adjoint values into the characterization of optimal control (54); 3. Repeat step 2 if a stopping criterion is not satisfied.