Optimal control of a dengue model with cross-immunity

Mathematical modelling of a dengue epidemic with two serotypes including a temporary cross-immunity yields a nonlinear system consisting of ordinary diﬀerential equations (ODEs). We investigate an optimal control problem, where the integral of the infected humans is minimised within a time interval. The controls represent human actions to decrease the number of mosquitos in the model. An integral constraint is added, which takes a limitation on the sum of the human actions into account. On the one hand, we derive and apply a direct approach to solve the optimal control problem. Therein, a discretisation of the controls is constructed using spline interpolation in time. Consequently, a ﬁnite-dimensional constrained minimisation problem can be solved. On the other hand, we employ an indirect approach, where necessary conditions for an optimal solution are considered. This technique yields a multipoint boundary value problem of a larger system of ODEs including adjoint equations. We present results of numerical computations, where the two methods are compared.


Introduction
Computational epidemiology helps to better understand the spread of diseases or the impact of a public health intervention, see [12], for example.Often compartment models are examined including populations of susceptible (S), infected (I) and recovered (R).Typically, the derived nonlinear dynamical systems consist of ordinary differential equations (ODEs) or delay differential equations (DDEs), see [1].
Concerning dengue fever, four serotypes of the virus are distinguished.An infection by some dengue serotype generates permanent immunity to it, but just temporary crossimmunity to other serotypes [5].SIR models have been designed for this class of problems, see [3,6].Furthermore, optimal control problems of dengue SIR models were examined concerning optimal vaccination strategies: with cross-immunity in [4] and without crossimmunity in [7].
We consider a dengue fever model used in [10], where two serotypes of the virus are included.The model consists of an initial value problem (IVP) of nonlinear ODEs describing the dynamics of ten human compartments and four mosquito compartments.In [10], an optimal control problem is investigated, which aims at decreasing the number of infected humans by actions to reduce the number of mosquitos.The objective function includes a combination of the number of infected humans and the costs of the used actions.Alternatively, we define an optimal control problem with an objective function depending only on the number of infected humans, whereas an additional constraint in integral form describes a strict limit on the total sum of the actions.
In general, an optimal control problem of ODEs can be solved numerically using a direct approach or an indirect approach, see [8].In the direct approach, a discretisation yields a constrained or unconstrained optimisation problem defined on a finite-dimensional domain.Thus numerical methods for classical optimisation problems generate approximate solutions.In the indirect approach, the local minimum principle typically implies a boundary value problem (BVP) of ODEs including adjoint equations.Consequently, numerical solvers for BVPs are applied.
We study both a direct method and an indirect method to solve the optimal control problem of the dengue fever model, where an IVP of the ODE model is considered.The controls are discretised using a finite number of time points in the direct approach.We construct a continuous approximation of the controls using a monotonicity-preserving Hermite spline interpolation.The Hamiltonian function associated to the system of ODEs is employed in the indirect approach.Necessary conditions for an optimal solution generate additional ODEs for Lagrange multipliers.Hence we obtain a (multipoint) BVP of ODEs with dimension twice as large as the original ODE model.Numerical results are presented for both approaches.We compare the properties and the performance of the numerical methods.

Problem definition
We consider a dengue fever model examined in [10], where two serotypes of the virus are included.The nonlinear dynamical system consists of 14 ODEs.There are 10 human populations and 4 mosquito populations, which are illustrated in Table 1.Let N h be the constant total number of humans.The dynamics of the human populations is modelled by and the dynamics of the mosquito populations reads as ( Initial values are predetermined at time t = 0.The right-hand side of the first equation in (2) has a term describing a logistic growth of the mosquito larvae.An optimal control problem of an ODE model including logistic growth was also investigated in [11], for example.In [10], the (non-negative) equilibria of the ODE system (1), (2) were determined in the case of constant parameters.There are up to two disease-free equilibria and up to two endemic boundary equilibria.The spread of the dengue disease is reduced by decreasing the number of mosquitos.The constants c a and c m represent actions to exterminate the larvae and the adult mosquitos, respectively.The capacity of the aquatic habitat is K m = k m N h .This capacity is decreased to αk m N h with 0 < α ≤ 1 by human actions.Now we consider these three possibilites as time-dependent controls c a (t), c m (t), α(t).In the following, we substitute the control α by the control c k satisfying Consequently, increasing c k causes a reduction of α.The range α ∈ (0, 1] corresponds to the range c k ∈ [0, ∞).The controls c a , c m and c k are restricted to c (t) ∈ [0, c ,max ] for t ∈ [0, T] and ∈ {a, m, k}.
( 3 ) In the system of ODEs (1), (2), the functions of the right-hand side are quasi-positive.It follows that all 14 populations remain non-negative for all t ≥ 0 provided that the initial values are non-negative, see [16, p. 84].This property is also valid in the case of our timedependent controls (even if the controls were negative).
As quantity of interest, we observe the total number of infected humans for t ∈ [0, T].It holds that I h (t) ≥ 0, since each population is non-negative.We want to minimise the number of humans infected within a time interval.Thus the integral of ( 4) is considered, i.e., We define an objective function by J = J(T).We can add the equation J = I h to the system of ODEs (1), (2).Hence a numerical integration of the enlarged system also computes the value J.
In [10], the objective function is minimised including the constant weights γ i , γ m , γ a , γ α > 0. The functional (6) takes the costs for the control of the mosquito populations into account.The term 1α(t) is used in (6), because some effort is required to make α(t) small, whereas increasing c a (t) or c m (t) produces expenses.However, there is no restriction on the magnitude of the controls.
A strongly increasing number of infected humans causes high control functions, and thus expensive costs have to be paid.Alternatively, we minimise the objective function J = J(T), see ( 5), subject to limitations on the control.We apply the constraint On the one hand, the constants γ m , γ a , γ k > 0 describe the price per unit of the different controls.On the other hand, the constant > 0 represents the available total funds or money for fighting against the disease of dengue fever.
In our numerical computations, we applied a scaled system of the ODEs (1), (2), where each equation is divided by the number N h .For example, the first equation of (1) changes into A scaling of the time variable is also used in the indirect approach.These transformations represent a nondimensionalisation of the system.Such nondimensionalisations could be applied to determine the relative importance of the terms in the right-hand side of the ODEs.However, this analysis is not within the scope of this article.
We suppose non-negative initial values, where the sum of all human populations is equal to one.It follows that the sum of all human populations is equal to one for all times.We apply the initial values whereas the other initial values are zero.Consequently, 99.99% of the humans are susceptibles at the beginning.The choice of the initial values concerning the mosquito populations follows [10].
We solve the IVP of the ODEs (1), ( 2) until a final time of T = 365 (days) in the case without control (c a , c m , c k ≡ 0).Table 2 lists the used parameter values.An implicit time integration using the trapezoidal rule yields a numerical solution of this IVP.Figure 1 illustrates the computed populations.Concerning the human populations, the plot includes the susceptible S h , the recovered R h , the infected I h from (4) and the sum of all others O h .Nearly all humans are infected by at least one of the two dengue types within the time interval [0, 100].In contrast, the mosquito populations A m and S m tend to stable nonzero equilibria.Alternatively, the aim of the optimal control consists in a reduction of the mosquito populations to decrease the total number of infected humans.

Direct approach
Now we derive a (so-called) direct approach to solve the optimal control problem of Sect. 2 numerically.

Discretisation
In the time interval [0, T], we arrange the grid points t i = (i -1) t for i = 1, . . ., M using a step size t = T M-1 .The associated discretisation of the controls reads as which implies 3M degrees of freedom.Let 20 ≤ M ≤ 40, for example.We construct continuous approximations c (t) for t ∈ [0, T] by an interpolation of the discrete values (8).
It turns out that a piecewise cubic Hermite interpolation including a preservation of monotonicity yields good results in this direct approach.In contrast, linear interpolation, which is monotonicity-preserving, or traditional cubic spline interpolation, which is not monotonicity-preserving, produce worse results.The continuous approximations c (t) = c (t; c ,1 , . . ., c ,M ) are inserted into • the right-hand side of the ODEs (2), and • the constraint (7) in integral form.The objective function J = J(T), see (5), also depends on the controls, since the controls affect the number of infected humans.Let x = (c a,1 , . . ., c a,M , c k,1 , . . ., c k,M , c m,1 , . . ., c m,M ) be the vector of unknowns.The discretisation yields a minimisation problem including an objective function J : R 3M → R on a finite-dimensional domain.Each evaluation of the function in (9) requires the solution of an IVP of the ODEs (1), (2).The box constraints (3) imply bounds on the discrete controls 0 ≤ c ,i ≤ c ,max for ∈ {a, k, m} and each i. (10) The integral constraint (7) changes into with a function G : R 3M → R. The constraint ( 11) can be evaluated with high precision using a relatively low computational effort.Now the minimisation problem (9) together with the constraints ( 10), ( 11) is solved by appropriate numerical methods, see [9], for example.Traditional numerical techniques are sequential quadratic programming (SQP), active set methods, interior point methods, and others.

Numerical solution of minimisation problem
The software package Matlab [14] was used for the numerical computations presented in the following.We consider the time interval [0, T] = [0, 365] in the optimal control problem.In the constraint (7) or (11), the weights are always set to γ m = γ a = γ k = 1, whereas the three cases ∈ {20, 30, 40} are studied.Furthermore, the bounds (10) are imposed with c ,max = 0.1 for ∈ {a, k, m}.We choose M = 30 in the discretisation (8).Consequently, the nonlinear minimisation problem includes 90 degrees of freedom.The constrained minimisation problem with objective function ( 9) is solved numerically using the Matlab built-in function fmincon, where an SQP algorithm performs the iteration.As starting values, we apply c (0)  ,i = 0.001 for ∈ {a, k, m} and all i = 1, . . ., M, (12) which results in a better performance than starting values identical to zero.
In each iteration step of the minimisation, IVPs of the ODEs (1), (2) have to be solved, which represent mildly stiff problems.Yet implicit integration methods show failures and sometimes produce bad approximations, where nonlinear systems are solved by Newton iterations.Thus we use the Matlab built-in function ode45, which executes an explicit Runge-Kutta method of order 4(5) including a step size selection based on a local error control.More details on numerical methods for solving IVPs of ODEs can be found in [17,18], for example.The limitations on the actions against the mosquitos are substantial in the case of = 20.
Figure 2 (left) shows the controls determined by the minimisation method.We observe that the controls c a and c k are nearly zero.Instead the control c m is significantly employed.
The ODE solutions using these controls are demonstrated by Fig. 3.In the case of no controls, all humans get infected within the time period.Now about 75% of the humans come down with dengue.
Alternatively, a plenty of actions against the mosquitos is possible in the case of = 40.3 itemises some statistics of the SQP iteration in the cases ∈ {20, 30, 40}.The starting value of the objective function ( 9) is always J = 5.0994.As expected, larger parameters yield lower final values of J.The final value of constraint (11) is shown, which is close to zero.Thus a maximum amount of actions is invoked in each case.

Indirect approach
We employ a (so-called) indirect approach based on necessary conditions for an optimal solution, see [2,15].

Boundary value problem
Representing the 14 populations of the dengue model by a vector-valued function y : [0, T] → R 14 , the ODEs (1), ( 2) can be written as Since the constraint ( 7) is linear in the controls c m , c a , c k , the augmented Hamiltonian represents a linear function, too: where H 0 = λ f 0 and H = λ f + νγ for ∈ {a, m, k} with some constant value ν and Lagrange multipliers λ : [0, T] → R 14 , including also the weights γ from (7).The necessary conditions for an optimal solution imply the system of ODEs together with the end conditions Finally, the optimal controls ĉa , ĉm and ĉk satisfy the maximum principle of Pontrjagin: for all admissible control values c a , c m , c k and for almost every t.According to the linear structure (13) and depending on y(t), λ(t), ν the control laws might be of bang-bang type again for ∈ {a, m, k} (16) or they might be singular, if H (y(t), λ(t), ν) = 0 for all t on some non-empty interval [τ 1 , τ 2 ].
The functions H are switching functions and a mapping t → (c a (t), c m (t), c k (t)) is called switching structure.Under the assumption that for some t ∈ [0, T] the mapping is known, namely c s is singular at t for all s ∈ S ⊆ {a, m, k}, the appropriate non-singular controls for all s / ∈ S could be determined subject to (16).In contrast, for all s ∈ S the controls c s satisfy the linear system of differential equations Due to the complexity of the underlying ODE system and the number of different control laws, the formulas for the singular cases are quite long and therefore not given here.The corresponding computations had to be done by the computer algebra system Maple [13].
Obviously, an appropriate switching structure is a priori unknown in our problem, as the switching functions H depend on the unknown values of y, λ and ν.

Numerical approach
The switching structure has to be guessed for a numerical solution.For this purpose, the time horizon is split up into several subintervals, such that within each subinterval the type of the controls is fixed throughout the subinterval.Hereby, the length of a subinterval becomes a free parameter of the solution process.In this way, we obtain a multipoint BVP, based on the ODEs (1), ( 2) and ( 14) together with the initial conditions for the state variables y and the end conditions (15) for the Lagrange multipliers λ.Since the degrees of freedom are the functions y, λ and the constant ν, the contraint ( 7) is added to our system of equations.Provided that the guess of the switching structure is convenient for the initial values of y and the physical parameters of the problem (Table 2), the BVP can be solved by some numerical method.We used the built-in routine bvp4c of Matlab [14], which implements a collocation method.
In the present problem, the switching structure depends significantly on the physical parameters (Table 2), the box-constraints (3), the weights and the upper limit in (7).The choice of suitable switching structures will be demonstrated subsequently for the value c m,max = 0.1 and different values of .We start by simplifying c a,max (t) = 0 and c k,max (t) = 0, for all t ∈ [0, T] ( and use = c m,max T = 36.5.In this situation, it is obvious that is optimal, since ( 7) is satisfied and the mosquito reduction is maximal on condition (3).The solution of the corresponding two-point BVP could be found and it yields the nega- tive switching function H m depicted in Fig. 5 (left).Hence the solution satisfies the condition (16).
For smaller values of it is reasonable to look for the switching structure including a free switching time τ .The switching functions of the solutions of the corresponding three-point BVP are shown for three values of in Fig. 5 (right) in the relevant time interval [300, T].The blue line is the same as the solution for (19), i.e., τ = T.The green line satisfies the necessary conditions of first order, in particular (16).The black line represents a solution of the three-point BVP, but it does not satisfy ( 16) within some intervals in some vicinity of the switching point τ = 340.Therefore another switching structure had to be found for this case.
Based on some experiments, a new switching structure for values ∈ [29.0, 34.0] is established that includes a singular phase: For = 34.0 the numerical solution satisfies all first-order conditions with a tiny singular phase (τ 2τ 1 ≈ 3.4).This interval enlarges with decreasing values of until = 29.0,see Fig. 6, where ( 16) is still valid, but the given box constraints (3) are violated on some subinterval of the singular phase.
Introducing another phase of maximum control inside the singular phase yields the next switching structure, which is appropriate within the parameter range ∈ [17.0, 29.0]:An initial phase of minimum control is added to the switching structure, to get Figure 8 (right) shows the results for the corresponding multipoint BVP in the initial phase with respect to the switching function H m .There is a very short minimum control phase (τ 1 < 1) depending on the value of the parameter .Figure 9 illustrates the solutions for the total time interval.The solution of the case = 15.5 (black line) indicates that again the control is not feasible close to t ≈ 40.An appropriate switching structure for this value could not be discovered.
A look at the switching functions of the neglected controls c a and c k reveals that, except for values of close to 36.5, these functions are positive throughout the entire time horizon.Therefore, in most cases investigated, the simplification ( 18) is not a genuine restriction in comparison to the original formulation.

Comparison of methods
Now we compare the performance of the two presented methods: the direct approach of Sect. 3 and the indirect approach of Sect.4, when applied to solve the optimal control problem under investigation.
In the direct approach, the minimisation problem (9) together with the constraints (10), ( 11) is solved separately for parameters = j and j = 1, 2, . . ., 40.We apply the numerical method as described in Sect.3. The starting values of the controls are always set to (12) in the iterations.The SQP iteration method using the scheme fmincon converges in each case to an approximation of a (local) minimum.Figure 10 illustrates the data of the objective function (9) for the approximate solutions.The value = 0 represents the special case of controls identical to zero, as examined at the end of Sect. 2.
In the indirect approach, the solution of the multipoint BVPs (1), ( 2),( 14) together with ( 16) and ( 17), respectively, is computed for parameter values ∈ [15.5, 36.5].The resulting data of the objective function is also shown in Fig. 10.Only the control c m is applied in these problems.We use a homotopy method (also called continuation method), see [18], with respect to the parameter .For the initial homotopy value ( = 36.5)the initial solution guess required by the collocation method bvp4c was chosen to be constant in all state variables, controls and adjoint variables.The number of collocation points is determined automatically by the scheme bvp4 and varies in our computations from about 400 to 900, depending on the size of the homotopy step and the value of the homotopy parameter.The initial solution guesses for the succeeding values of are taken from the preceding steps.The objective function is constant for ≥ 36.5, because the maximum amount of c m is employed due to the box constraint (3).Thus the controls c a and c k would have to be included to decrease the objective function further.The homotopy method stops at = 15.5, because the numerical solution of the multipoint BVPs fails for lower values of the parameter.In total, 252 BVPs are solved by the homotopy method, where a non-constant step size is used for the parameter.
Concerning the final values of the objective function, the results from the indirect approach are always better or equal than the results from the direct approach.However, the differences are relatively small between the two methods.Moreover, the computational effort is much larger in the indirect approach.Furthermore, the direct approach does not require the usage of a homotopy method.The direct technique yields good approximations for all parameter values, whereas the indirect technique fails for small parameter values.Hence the direct approach exhibits better convergence properties in this ODE problem.
Finally, we comment on the behaviour of the minimised objective function in dependence on the parameter , see Fig. 10.Starting from small values of , the objective function decreases nearly linearly.There is a critical range around ≈ 30, where the objective function decreases rapidly.Now the available actions, which are described by the controls, are sufficiently extensive such that the disease is nearly exterminated.

Figure 2 Figure 3
Figure 2 Controls computed by direct method in the two cases = 20 (left) and = 40 (right)

Figure 2 (
Figure 2 (right) depicts the controls computed by the minimisation method.Still predominantly the control c m is activated, which reaches the upper bound most of the time.Yet the controls c a and c k are also mobilised at the beginning of the time interval.The ODE solutions associated to these controls are shown in Fig. 4. Now nearly no humans incur an infection during the time period.Although the number of mosquitos is still significantly larger than zero, the number of infected mosquitos approaches almost zero.Moreover, we compare the performance of the minimisation methods for different parameter values.Table3itemises some statistics of the SQP iteration in the cases ∈ {20, 30, 40}.The starting value of the objective function (9) is always J = 5.0994.As ex-

Figure 4
Figure 4 Human compartments (left) and mosquito compartments (right) from ODE model using controls in the case of parameter = 40

Figure 10
Figure 10Objective function for different parameter obtained in optimal control problems solved by direct approach or indirect approach (left: linear scale, right: semi-logarithmic scale)

Table 2
Parameter values in dengue fever model 3 Figure 1 Human compartments (left) and mosquito compartments (right) from ODE model in the case of no control

Table 3
Performance of iterative minimisation method for different parameter values of