Model-Based Monitoring of Dengue Spreading

The problem of designing a system for monitoring the spread of dengue in a mixed population of humans and mosquitos is addressed. For this purpose a model–based observer approach is employed based on a model that describes the actual number of infected humans (with constant population size) and mosquitos (with time–varying population size) and which has been previously validated with data from different regions of the state of Morelos in Mexico using the number of new reported cases as available measurement. This model has time–varying reproduction parameters for the mosquitos and thus enables to account for different climatic conditions in the course of the year (in particular the rain and dry periods) that affect the mosquito population size. An analysis of the structural observability, based on the preliminary assumption of a continuous measurement of the accumulative number of reported cases of human infections is carried out to motivate a reduced order model that is structurally observable. For this reduced order model a continuous–discrete extended Kalman Filter is designed as basis for the monitoring and prediction scheme. The approach is validated using real data from Cuernavaca, Morelos, Mexico, showing a high potential for future developments towards automated monitoring schemes.


I. INTRODUCTION
Dengue fever as an infectious disease caused by the dengue virus and transmitted to humans by infected mosquitoes of the genus Aedes aegypti that affects millions of people in tropical and subtropical regions. According to estimates by the World Health Organization (WHO) [1] about 390 millions people get infected every year out of which only 96 million receive medical treatment. Due to the huge population affected by the virus, the WHO has classified this disease an epidemic. For the effected countries, mainly in South Asia and South America it additionally implies a considerable economical burden [2], [3].
Currently, there are several approaches for predicting Dengue spreading in endemic regions of the world. A rather complete literature review was performed by Sylvestre et al. [4], with most of them being based on machine learning The associate editor coordinating the review of this manuscript and approving it for publication was Angel F. García-Fernández . approaches (see also [5], [6], [7]). An interesting survey of different mathematical models for dengue's spreading prediction is given in [8]. From a population dynamics perspectives different approaches have been proposed for the modeling going from the causal [9], [10], [11], [12], [13], [14], [15] to the statistical modeling [16], [17].
In several recent studies [18], [19] the problem of dengue outbreak estimation and forecast was addressed using some types of the vector-borne disease model SEIR-SEI with compartments Susceptible-Exposed-Infectious-Recovered (for host) and Susceptible-Exposed-Infectious (for vector) in combination with ensemble Kalman filtering (EnKF) for data assimilation. These studies highlight the great potential of combining compartmental models with Kalman Filters for achieving high-accuracy now-casts of non-directly measured states and on this basis being able to provide better forecasts for the upcomming weeks and months.
In particular, understanding of the mechanisms involved in the spreading not only in the human population but also in the FIGURE 1. Classical SIR state transition diagram for the human population with susceptible (S), infectious (I) and recovered (R) states and associated transition rates τ i →j . mosquito population can provide means for better designing monitoring and forecasting (i.e., early warning) systems for prevention, containment and possible on-time control of the spread of the disease (see, e.g., [13], [20]).
The main purpose of the present study is to further analyze and employ the model proposed in [15] with respect to its structural observability properties and the design of a monitoring and prediction scheme based on a nonlinear state observer. Based on the observability analysis a reduced order structurally observable model is determined. For this reduced model a continuous-discrete Extended Kalman Filter [21] is designed and its performance is established through numerical simulations with a parameter set that has been validated in [15] for the city of Cuernavaca, Morelos, Mexico. The potential use of this state estimation technique for prediction of the future spreading is shown using actual data from Cuernavaca from the year 2012.
The paper is organized as follows. In Section 2 the monitoring problem is discussed and the mathematical model presented. The analysis of the structural observability is carried out in Section 3. The design of the continuous-discrete Extended Kalman Filter is presented in Section 4. In Section 5 the validation and testing by numerical simulations is presented. Final conclusions and the future outlook are presented in Section 6.

II. MATHEMATICAL MODELING
Consider a human population with N ∈ N inhabitants, each of which can be in either of the states S (susceptible), I (infectious), R (recovered), according to the state transition diagram shown in Figure 1.
In this state automaton the state transition rates τ S→I , τ I→R are both positive. Given that dengue is not transmitted between humans and people get infected only due to the bite of an infected mosquito and are assumed to recover after a certain period (implying a constant population) and become later on again susceptible, given that there is no immunity against the dengue virus. Given that infections are triggered from outside the human population a possible approach in the modeling consists in considering the mosquito population an exogenous time-varying disturbance input affecting the dynamics of the states S, I, R and thus the reports of total accumulated infection numbers. This approach is schematically depicted in Fig. 2.
This approach could be sufficient for the monitoring of the human population using, e.g., an adaptive observer.
An alternative approach, providing more physical insight into the underlying dynamics by representing further   important mechanisms in the spreading process consists in providing a model for the spreading within the mosquito population. For this purpose consider a time-varying population of mosquitoes which can be in either of the states U (healthy), V (incubating), W (infectious). In contrast to the human population, mosquitos do not live sufficiently long to recover, once they are infected, but can indeed pass the infection to their offspring. This yields the following state transition diagram shown in Figure 3 for the mosquito population with the positive state transition rates τ R→S , τ U →V , τ V→W . As their are clear interactions between the mosquito and human populations one obtains a feedback loop as shown in Fig. 4 for the overall spreading process dynamics.
Denoting in the following the number of susceptibles, infectious, and recovered humans by S, I , R, respectively, the number of susceptible, incubating and infectious mosquitos by U , V , W , respectively, and considering mass-action-like transition rates it holds that In the preceding equations N m denotes the number of mosquitos, α represents the probability that a bite from an infected mosquito renders a susceptible human infected, and δ the probability that when a susceptible mosquito bites an infected human the mosquito gets infected. The terms W /N m and I /N determine the probability that an arbitrarily chosen mosquito is infectious or an arbitrarily chosen human is infectious, respectively. Having these transition rates the VOLUME 10, 2022 total rate of change in S and U due to bites can then be computed by τ S→I S and τ U →V U , respectively. With the number of accumulated cases denoted as A, the dynamics of the spreading process between both populations can be modeled by the equation seṫ for t > 0 with initial values In the preceding equation set the parameter γ > 0 represents the rate with which recovered humans loose their immunity, ρ > 0 the recovery rate, κ > 0 the rate with which incubating mosquitos become infectious, r 0 , r 1 , r 2 , r 3 ∈ R are the net birth and decease rates. The parameter K represents the (constant) maximum mosquito population size. It is assumed that S, I , R, A, U , V , W : [0, ∞) → R + are continuously differentiable functions with R + denoting the positive real numbers. Given that the right-hand side of the differential equation system (2) is locally Lipschitz continuous, a unique solution exists that is continuously depending on the initial conditions.
Assuming a constant population of humans the total population numbers are given by Note that the first equation corresponds to an algebraic constraint that is naturally satisfied as long as S 0 + I 0 + R 0 = N . This is summarized in the following lemma. Lemma 1: [15] The total human population size N = S + I + R remains constant over time.
The next lemma shows an interesting property of the model (2) putting it in direct relation to classical population balances as, e.g., the Verhulst model [22].
Lemma 2: [15] For constant reproduction parameters r 1 = r 2 = r and r 3 = r − r 0 the total mosquito population size N m = U + V + W is determined by the logistic growtḣ For non-constant reproduction parameters r k , k = 0, 1, 2, 3 the mosquito population will be time-varying. Note that, as the reproduction rates will depend on the climatic conditions [23] it makes sense to consider time-varying rates. One simple way to take seasonal differences in reproduction rates into account, which mostly are present due to rain and dry seasons, the birth and decease rates of the mosquitos are considered piecewise constant according to for k = 0, 1, 2, 3. The time varying rates and associated non-constant mosquito population size with Verhulst-like dynamics represent important differences in comparison to other models that have been reported in the literature (see, e.g., [12] and references therein).
In compact form the model (2) is written aṡ

III. STRUCTURAL OBSERVABILITY ANALYSIS
In this section the structural observability [24], [25], [26] of the spreading model (2) is analyzed in the understanding that the structural observability is a necessary condition for the (local) complete observability and provides an important basis for design of an observer-based monitoring scheme. Definition 1: The system is structurally observable if it is completely observable for at least one state-parameter pair.
For the analysis of the structural observability first recall the definition of the structure graph (x) = (V , E(x)) of the system at a point x ∈ X , consisting of the vertex set V = {v 1 , . . . , v 7 } where v k is associated to the state x k , k = 1, . . . , 7, and the (directed) edge set ∂x j = 0}. Note that the edge set in principle depends on the value of x, given that ∂f i (x) ∂x j can be zero at some particular states. Structural observability is a generic property for a system that depends only on the intrinsic interconnection structure and not on the particular values of ∂f i (x) ∂x j ∈ R, i.e., it is actually independent of the particular parameter set and the particular value x where the partial derivatives are evaluated. It only depends on whether the influence of x j on x i is not always zero. Indeed, if it is not zero for at least one parameter vector and one state value x ∈ X , due to continuity reasons it will be nonzero for almost all values of x and parameters, up to those lying on a specific hyperplane of the state-parameter space which has measure zero (see, e.g., the discussion in [25]).
The characterization of possible graphs implying structural controllability and observability can be reviewed, e.g., in [24], [26] and typically implies the extension of the state graph discussed above by input and output vertexes. As in the present study no inputs are considered and the measured FIGURE 5. Structual connectivity graph for the 7-state bi-species model with measurement of the state A. The red dotted connection from S to R shows the lack of structural observability, given that from R one can not reach any other node without crossing I again, implying a loop.
output corresponds with the state A the essential information is already contained in the state graph. The condition for structural observability then basically requires that a unique information path exists from the measured nodes toward all other nodes without loops and intersections. For the case at hand this can be readily applied to the structure graph shown in Fig. 5. Based on this graphical representation the following result is obtained.
Proof: As can be seen in Fig. 5 the signal flow cannot pass through the complete graph toward the measured node A without passing twice through the node S, implying that no information path without loops exists.
As a particular outcome of this analysis one can see that it is possible to obtain a straight information path without loops from the output node A to all other nodes except R, as highlighted by different colors (and shapes) of the connections in Fig. 5. Recall from Lemma 1 that the value R can be obtained at any moment from the algebraic relation and one can reduce the model by this state without loosing any information. This leads to the reduced order model This model is identical to (2) with R in the first equation substituted by the relation (7). The structure graph of this model is shown in Fig. 6 and the following result is an immediate consequence of the above considerations. Lemma 4: The reduced order model (8) is structurally observable.
In compact form the reduced model (8) is written aṡ with the state vector x r (t) ∈ R 6 , where

IV. CONTINUOUS-DISCRETE OBSERVER DESIGN
Having the continuous-time model (8) (or equivalently (9)) that is structurally observable it would be straight forward to implement a continuous-time observation scheme. Anyway, as the reporting of new cases is carried out in discrete, periodic 1 time steps a continuous-discrete observation scheme is considered in the following. In this scheme the continuous-time model is used for predicting the state in between sampling times and the correction only takes place at the sampling moments. A well-established approach to address this kind of problem is the continuous-discrete Extended Kalman Filter (cd-ekf) [21]. Given that different stochastic sources are influencing the spreading process which are not modeled in (8) additional stochastic perturbation terms are included (so-called process noise) and measurement uncertainties are summarized in a measurement noise component. Assuming the process and measurement noise variable w and v, respectively, to be normal distributed with zero mean and covariances Q, R, respectively, this yields the stochastic model The basic idea consists in using (i) the deterministic nonlinear reduced model (8) for predicting the state between measurements, (ii) a linearization around the actual estimated state for the prediction of the associated estimation error covariance, and (iii) a correction at the discrete measurement instances t k , k ∈ N of the state and covariance estimates. This scheme is summarized in the following algorithm. a) Initialization: Set k = 0,x r (0), P(0). b) Prediction: Set the initial conditionsx r,p (t k ) =x r (t k ) with the estimatex r (t k ) of the reduced state x r (t k ) in (9), P p (t k ) = P(t k ) with the estimate of the covariance matrix, and solve the reduced differential equation systeṁ over the time horizon t ∈ (t k , t k+1 ] to obtain the predictionsx r,p (t) and P p (t). c) Correction: Update the predicted value with the actual measurement y(t k+1 ) = A(t k+1 ) + v in dependency of the predicted estimation error covariance P p (t k+1 ) as with the Kalman correction gain L k+1 . Repeat from b) with k = k + 1.
Generally speaking, the aim of the Kalman Filter is to minimize the estimation error covariance P(t) by appropriately choosing the correction gain L k (see, e.g., [21]).
At any time instant t ≥ 0 the estimate of the number of recovered humans is given bŷ It should be mentioned that besides the cd-ekf proposed here, e.g., unscented [27] or ensemble-based [28] Kalman Filters, particle filters [29], or nonlinear observers could be employed for the task of providing an estimatex r .

V. VALIDATION
In this section the performance of the proposed approach is tested for the problem of monitoring the spread of dengue in the state of Morelos, Mexico within the city of Cuernavaca for which monthly reported new cases have been used to determine the increasing number of accumulated cases A. In [15] a parameter set for the spreading during the year 2012 has been presented and is given as follows: , For the simulation the following initial condition for the system (x 0 ), Kalman Filter (x 0 ), estimation error covariance (P 0 ), variance of measurement noise (R), and covariance matrix for the process noise (Q) have been used: x(0) = 369996 4 0 4 10 4 6.19 109 x r (0) = 369930 54 14 20 1500 9.285 P 0 = I , R = 1 Q = diag 10 3 10 2 10 2 10 3 10 3 10 3 , The associated time evolution of A andÂ are shown in Fig. 7, and the resulting state evolution for S, I , R, U , V , W and their estimatesŜ,Î ,R,Û ,V ,Ŵ are shown in Fig. 8. It can be seen in Fig. 7 that the value of A is quickly recovered at the measurement time instances t k , k = 2, . . . , 12, with an increasing deviation between measurement time instances due to wrong estimates in the remaining states, but that after about 4 months (i.e., 3 measurements) the estimated valueÂ almost exactly follows the simulated value A. For the states, one can see in Fig. 8 that the states associated to the human population (i.e., S, I ) are rather well recovered, while the ones for the mosquitos are only qualitatively reconstructed. The reason for this can be manifold and will need to be investigated with more detail in future studies. Given that the aim of the proposed monitoring scheme is to recover the information about the human infection process and does not focus on the estimation of the mosquito population, with the latter only being included to cover the main mechanisms involved to retain the underlying feedback structure, the result seems reasonably well working to be implemented as a monitoring mechanism.  . Example of prediction from the actual estimate at time t = 5 with data (stars), uncorrected prediction starting from t = 4 (red dotted line), corrected prediction starting from t = 5 (yellow dash-dotted line) and the associated ±3σ A variance band (shaded region starting from t = 5). Time instances t < 5 belong to the past and time instances t > 5 to the future.
A particular advantage of the corrections introduced at measurement time instances is that they improve the prediction capabilities of the model. To highlight this in Fig. 9 a comparison is provided between predictions based on the initial value at time t = 5 for the interval [5,12] without correction (red dotted line) and with the correction using the measurement at t = 5 (yellow dash-dotted line). The stars indicate the actual reported cases numbers for Cuernavaca during 2012 (starting from February). It can be clearly seen that an advantage is obtained by using the automated correction for improving the estimate at t = 5 and using the corrected estimate for future predictions. In the figure there is also the ±3σ A interval shown around the prediction value, with σ A being the predicted estimation error variance for the state A, comming from the approximated covariance evolution dynamics (12b) (without further corrections after t = 5). This highlights the potential of using this approach including automated correction with new reported cases for the monitoring and prediction of the spreading of dengue.

VI. CONCLUSION
The problem of monitoring the spread of dengue in a mixed human-mosquito population has been discussed. A monitoring scheme has been designed that is based on a recently proposed validated mathematical model. Motivated by a structural observability analysis a reduced order model has been used for the design of a continuous-discrete Extended Kalman Filter. The proposed estimation scheme has been tested with the validated model for a case study for actual reported cases in the city of Cuernavaca, Morelos, Mexico. The potential use of the proposed scheme for monitoring and prediction based of improved initial estimates provided by the Kalman Filter has been discussed and illustrated.
Future work will address further validation studies based on different data sets, as well as further analysis of the underlying model in terms of detectability, non-local observability and the design of robust observers. Given the intrinsic dependencies of the reproduction of mosquitos on climatic conditions it is planned to combine these approaches with suitable machine-learning techniques for automated adaptation of the model parameters.