Electronic Journal of Qualitative Theory of Differential Equations Functional Continuous Runge–kutta–nyström Methods

Numerical methods for solving retarded functional differential equations of the second order with right-hand side independent of the function derivative are considered. The approach used by E. Nyström for second-order ordinary differential equations with the mentioned property is applied for construction of effective functional continuous methods. Order conditions are formulated, and example methods are constructed. They have fewer stages than Runge–Kutta type methods of the same order. Application of the constructed methods to test problems confirms their declared orders of convergence.


Introduction
Deriving mathematical models of epidemic propagation was primarily and originally motivated by the need of controlling the spread of a disease.This demand has led to several modeling approaches starting from compartmental models to network-based models with different complexity, such as mean-field, pairwise, compact pairwise, degree based or individual based models [1,10,13,15].All of these models are systems of non-linear differential equations, placing the problem of epidemic control into control theory [8,19], a well-developed Corresponding author.Emails: bodoagi@cs.elte.hu(Á.Bodó), simonp@cs.elte.hu(P.L. Simon) field of mathematics motivated mainly by engineering problems.Similarly to other questions, the control of epidemics can be formalised in the language of optimal control.The target is to drive the number of infected individuals to zero (or at least decreasing it as much as possible).This could be simply achieved by separating all the individuals from each other or vaccinating or treating all of them.This intervention obviously generates a cost that can be quantified by a cost function, the value of which is minimised.
Most models use vaccination or treatment as control measures, hence the cost function incorporates the cost of the vaccine or medical treatment [3,9,12].Since the literature of epidemic control is extremely rich, we mention here briefly only those results that are significantly relevant from the point of view of our study.The optimal time dependent vaccination was investigated in a susceptible-infected-recovered (SIR) model under minimising a cost function that measures the cumulative amount of infected and vaccinated people [16].Clancy and Piunovsky [3] computed optimal control based on isolation in a variant of the classical SIR model with nonlinear infection rate function.Hansen and Day [9] considered optimal control in the presence of limited resources using isolation, vaccination and mixed control strategies for the SIR dynamics.The application of optimal control theory to SIR propagation on networks in a real world situation is presented in [17].The effect of vaccination in an SIRS model on heterogeneous networks is studied in [2].These models are typically based on compartments of individuals in different states of the disease, in different age or place of working and living.Quarantine and contact tracing offer another control measure that uses information about the relation among the individuals in a certain sense [11,14].
An alternative is to incorporate the contact structure of the individuals and use the deletion or creation of links between them as control measure [18].This has been carried out for the pairwise SIS (susceptible-infected-susceptible) model in [18], where nonlinear model predictive control was applied to drive the number of infected nodes in the network to zero while keeping the network well connected.In that paper, the SI edges (connecting a susceptible and an infected node) are cut to eradicate disease and SS edges are created to maintain a given average degree in the network.
In this paper, we make a further step towards network-based epidemic control.The differential equation models of epidemic propagation on networks are approximations of the real process since all of them contain closure relations.Thus it is natural to ask how the original process can be controlled.The full stochastic model has an enormously large state space (for SIS epidemic on a network with N nodes, the state space has 2 N elements), hence individual based stochastic simulation is used to follow the process.Our goal in this paper is to control the stochastic simulation itself.We will study to what extent the control predicted by a network-based reduced ODE model results in good control of the full stochastic model.Studies in this direction already exist, and the first signs are positive.Namely, control computed from mean-field models seem to translate well, at least for some cases, for the control of the stochastic simulation counterpart [4,20].Our approach is motivated also by the fact that the cost of computing control from ODEs is much smaller than that of working out control from stochastic models.The main idea of our approach can be summarised briefly as follows.The process is controlled by cutting SI and creating SS edges in the stochastic simulation.The number of edges to be created or to be deleted is determined from the pairwise ODE model by updating its initial condition at the time instants k∆t, with an appropriately chosen time step ∆t, by using the values obtained from the stochastic simulation.
The results of our study are the following.On the one hand, we found that computing the values of the control parameters from the ODE approximation can be used to control the stochastic simulation.On the other hand, the success of the control strongly depends on the epidemic parameters, the infection and recovery rate, and on the control parameters, the length of the time horizon, the time step ∆t, of updating the control parameter values, and the parameters of the optimisation process in the model predictive control.We will show numerical evidence that the stochastic simulation can be controlled by using the control signals obtained from the ODE approximation, if the time step ∆t is small enough.The effect of the control bounds on the controllability will also be studied in detail.
The paper is structured as follows.In Section 2 our stochastic model for SIS epidemic propagation on an adaptive network and its pairwise ODE approximation are introduced.The detailed description of the control method is given in Section 3 both for the stochastic simulation and for the pairwise ODE.The results about the controllability of the simulation are presented in Section 4.

Model formulation
In this section the mathematical model of SIS epidemic propagation on an adaptive network with link-type dependent link activation and deletion is presented.To describe the model consider an undirected simple graph with N nodes, where each node can be either susceptible (S) or infected (I).A susceptible node can become infected when contacted with an infected one, and an infected one can recover and become susceptible again.In order to avoid infection, it is reasonable to assume that susceptible nodes break their links to infected ones and reconnect to a randomly chosen susceptible node, moreover susceptible nodes may cancel their links to other susceptible ones in order to avoid a possible infection if the neighbour becomes infected.Hence in our model it is possible to create new links and to delete existing connections leading to changes in the network structure.
In the standard adaptive SIS model, infection, recovery, link activation and link deletion are governed by independent Poisson processes.A susceptible node becomes infected at rate kτ, where k is the number of its infected neighbours and τ > 0 is a given constant, called infection rate.Each infected node recovers at rate γ, where γ > 0 is called the recovery rate.A non-existing link between two nodes of type A ∈ {S, I} and B ∈ {S, I} is activated at rate α AB .Similarly, an existing link between a node of type A ∈ {S, I} and a node of type B ∈ {S, I} is terminated at rate ω AB .The graph is assumed to be undirected, therefore we assume that α SI = α IS and ω SI = ω IS .
The mathematical model of the process is a system of ordinary differential equations, called master equations that consists of differential equations for the time dependent probabilities of each state.The state space of the adaptive SIS model on a networks with N nodes contains elements, because every node can be either susceptible or infected and every link can be either existing or non-existing.Therefore, even writing down the differential equations is not feasible for large N. To demonstrate the complexity of the problem let us take a look at the case N = 2. Let X AB denote the probability of the state, in which the first node is of type A and is connected to the second one, which is of type B. Similarly, the probability of the state where there is no link between the nodes is denoted by X AB .The master equations take the following form: where dot denotes differentiation with respect to time.
In order to overcome this difficulty, approximating non-linear differential equations have been derived for coarse-grained variables.This can be done at the individual level or at the population level.One of the latter is presented in the next section.

Pairwise model
A widely used population level approximation is the so-called pairwise model.This is formulated in terms of the number of nodes and edges in different states, hence it is suitable for incorporating the control action as creating and deleting certain types of links.The pairwise model for the expected values of the number of nodes and links takes the form (see Chapter 8 in [15]) (2.1) where N is the size of the population and n = n(t) is the actual mean degree of the network, i.e.

Network-based stochastic model
The stochastic simulation is based on the Gillespie algorithm ([6], [7]), which we have implemented in MATLAB.Let [0, T] be the time interval of the simulation.Moreover, let G = (V, E ) be the initial graph, where V and E represent the nodes and the edges of the graph, respectively.A node v ∈ V can be susceptible or infected and an edge e ∈ E can be of type SI, SS or I I. Let us assume that the parameters N = |V |, τ, γ, α SS , α SI , α I I , ω SS , ω SI and ω I I are all given.
The simulation is based on an accurate book keeping of all possible events, i.e. a node may become infected or susceptible, or an edge is created or removed.The most important step of the algorithm is to compute the rates of all possible events.These rates can be represented with a vector r of size where the first N elements correspond to the events related to the nodes (infection and recovery) and the remaining N(N − 1)/2 elements of r refer to the edges (creation and deletion).For instance, in case a node v is infected, i.e. it may recover, the rate corresponding to this node is γ.In case v is susceptible with k infected neighbours, its corresponding rate is kτ.An existing link AB ∈ {SS, SI, I I} can be removed, therefore the corresponding rate is ω AB .Furthermore a non-existing link AB ∈ {SS, SI, I I} can be created, therefore the rate belonging to it is α AB .Let R denote the sum of the rates, i.e.R = ∑ i r(i).The next step is to specify the time t of the next event.To determine this we choose a random number from an exponential distribution with parameter R. Next we determine which event will be performed.The event is chosen randomly but proportionally to its rate and may correspond to either a node or an edge.Once the event is chosen we update the graph G, i.e. the states of the nodes and the edges.After that using the new graph G the rates of the events are calculated, then a new time step t is determined and the next event is chosen randomly again.
Given that all existing and non-existing links need to be accounted for, the algorithm which includes the storage, update and referencing back and forth between rates and events becomes more complex.

Dynamic control of the adaptive SIS model
Now we turn to the main goal of the paper, to investigate the controllability of the adaptive SIS epidemic on a network by creating and deleting edges.Our aim is to eradicate the epidemic and keep the network well connected as well, i.e. to lead the system to the state when the expected number of infected nodes is [I](T) = I * and the average degree is n(T) = n * for some finite T > 0 and target values I * = 0 and n * = n(0).
In order to achieve our goal we change the rates of creation and deletion of SS links and deletion of SI links from time to time.This process models that healthy people terminate their connections to infected ones and try to find connections to other healthy individuals in order not to be disconnected socially.Thus in our control process we have and where u 1 and u 2 are our control parameters.Since we use time dependent control, u 1 and u 2 will be piecewise constant functions.We fix a step size ∆t > 0 determining how often we intervene and change the values of the control parameters.That is the control functions u 1 and u 2 will be constant in each time interval [(k − 1)∆t, k∆t), for k = 1, 2, . . ., K with K∆t = T. Furthermore we assume that the control parameters are bounded, namely there exist positive numbers M 1 and M 2 , such that We define the controllability of our system as follows.In practice, the case ε = 0, i.e. perfect controllability is not expected.
In the next subsections we show how the control functions u 1 and u 2 are determined for the pairwise model and for the stochastic simulation.

Time dependent control of the pairwise model
Using the parameter values given in (3.1)-(3.2),system (2.1)-(2.4)takes the following form: (3.7) Our aim now is to determine the values of the piecewise constant functions u 1 and u 2 in such a way that system (3.4)-(3.7)becomes ε-controllable.This has already been carried out in [18] by using nonlinear predictive control.Here we only briefly summarize the main steps of the method because these control signals will be used to control the stochastic simulation in the next subsection.The Nonlinear Model Predictive Control (NMPC), which is a variant of Model Predictive Control (MPC), is an optimization based method for the feedback control of nonlinear systems.It is widely applied to stabilization and tracking problems, see [8].
First, we discretize system (3.4)-(3.7) in order to apply the NMPC algorithm.In order to make the equations simpler use the following vectorial notations: Let us denote the output variables, i.e. the number of infected individuals and the mean degree, by Y = ([I], n).We use X i and Y i for the i-th coordinate of X and Y respectively, and x(k) = X(k∆t), y(k) = Y(k∆t), u(k) = U(k∆t).Let us denote the solution operator of system (3.4)-(3.7)by φ, that is the solution starting from the initial condition X(0) is given by X(t) = φ(t, X(0), U(t)).Now, the control is constant in a time interval of length ∆t, hence we will use the discretized solution function where x ∈ R 4 and u ∈ R 2 are given vectors.Furthermore, we introduce the discretized output function as the first coordinate of which is [I] and its second coordinate is the average degree.Following [18], we compute the control action u at time k∆t in the following way.First, we set a prediction horizon of length P steps and denote by u i (k + j) for i = 1, 2, j = 0, 1, . . ., P − 1 an arbitrary admissible future control action chosen at time k∆t.Given this control sequence, the solution at the end points of the time steps starting from x(k) = X(k∆t) is given as The output variables are simply given as y(k + j) = H(x(k + j)), for j = 1, 2, . . .P. Now comes the most important step, to compute the sequence of control values u(k), u(k + 1), . .., u(k + P − 1) in such a way that the output values y(k + 1), y(k + 2), . .., y(k + P) get as close to the target value (I * , n * ) as possible.This can be achieved by minimising the objective functional J : R 2P → R given as J (u(k), . . . ,u(k + P − 1)) where λ 1 , λ 2 , λ 3 , λ 4 are parameters (damping parameters) and ∆u i (k + j) = u i (k + j) − u i (k + j − 1) is the change of control value from the k-th time step to the next one.We solve the above nonlinear optimization problem using the so-called lsqnonlin optimization routine in MATLAB, yielding values for the control values u(k), u(k + 1), . .., u(k + P − 1).Then the control signal u(k) is used in the time interval [k∆t, (k + 1)∆t) and the whole computation is repeated at time (k + 1)∆t to obtain the control signal in the interval [(k + 1)∆t, (k + 2)∆t) and so on until the end of the time interval, T, is reached.
It was shown in [18] that choosing the control parameters appropriately, system (3.4)-(3.7) is ε-controllable for small values of ε, i.e.Y(T) will be sufficiently close to the target value (I * , n * ).In the next subsection, this result is applied to control the stochastic simulation.

Control of the network-based stochastic model
Now our aim is to control the network-based stochastic model.The question is how to choose the control parameters in order to control the simulation.The idea for controlling the pairwise model was to determine the output for an arbitrary admissible control sequence and then minimising the objective function to get the best control signal.In the case of the stochastic simulation, the difficulty is in computing the output for several different control sequences.Hence the idea is to use the control values obtained for the pairwise model in order to control the simulation.

Control of the stochastic model by using method (B)
Now we present successful control results obtained by applying the second method.We recall that in this case the control signal in the time interval [(k − 1)∆t, k∆t) is computed from the pairwise model but the initial condition at time (k − 1)∆t is taken from the stochastic simulation.
A successful control process is shown in Figure 4.2, where the time dependence of the prevalence, average degree and the control signals are shown.One can see that the link deletion rate, u 1 is at the maximum value, M 1 in the first part of the process.The creation and deletion of SS links, u 2 , shows an interesting time dependence.First, it is negative, that is SS links are deleted (in order to prevent strong propagation), then SS links are created in order to achieve the target value of the average degree.We note that this seems to be the result of intentional planning, however, this is simply given by minimising the objective function.This numerical experiment is repeated with a longer time step, ∆t.As it is expected, the control is not effective when the time step is not small enough, that is when the initial condition of the pairwise model is rarely updated from the stochastic simulation.In Figure 4.3 the time evolution of prevalence, network connectivity and control signals are shown for the same parameter combination as in Figure 4.2, but with ten times bigger time steps.One can see that the average degree does not reach its target value.Now we will investigate the effect of the control parameter M 1 with fixed values of τ and M 2 .In [18] it was shown that for a fixed value of τ and M 2 there is a critical value M c 1 such that if M 1 is below this critical value, then the control is ineffective, while for M 1 values above M c 1 , the control is effective for the pairwise model.It turned out that a similar threshold value of M 1 exists when the stochastic simulation is controlled, but in this case the threshold depends also on the length of the time step, ∆t.Fixing a value of M 2 and ∆t, we determined the threshold of M 1 for several values of τ.The results are shown in Figure 4.4.For a given value of τ, the curves yield the threshold of M 1 , that is the control is effective when M 1 is above the curve.
We observed that for several parameter combinations, also for those shown in the figure, the critical value of M 1 in the pairwise model is larger than that for the stochastic simulation, which means that the control of the stochastic model is "cheaper" than the control of the pairwise model.
[10]em (2.1)-(2.4) is not self-contained in this form, since there are no equations for the number of triples.Instead of deriving differential equations for those variables, an approximating algebraic relation is used typically to express the number of triples in terms of the pairs and singles.The most frequently used and widely accepted approximation is the following[10]: )where[I](t) and [S](t) are the expected values of the number of infected and susceptible nodes at time t, ([S](t) + [I](t) = N holds for all t), [SI](t), [I I](t), [SS](t) denote the expected number of SI, I I and SS edges and [SSI](t), [ISI](t) denote the expected number of SSI, ISI triples at time t.

Table 4 . 1 :
Difference between the prevalence and average degree and their target for different values of τ, M 1 and M 2 for N