Measles dynamics on network models with optimal control strategies

To investigate the influences of heterogeneity and waning immunity on measles transmission, we formulate a network model with periodic transmission rate, and theoretically examine the threshold dynamics. We numerically find that the waning of immunity can lead to an increase in the basic reproduction number R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0}$\end{document} and the density of infected individuals. Moreover, there exists a critical level for average degree above which R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{0}$\end{document} increases quicker in the scale-free network than in the random network. To design the effective control strategies for the subpopulations with different activities, we examine the optimal control problem of the heterogeneous model. Numerical studies suggest us no matter what the network is, we should implement control measures as soon as possible once the outbreak takes off, and particularly, the subpopulation with high connectivity should require high intensity of interventions. However, with delayed initiation of controls, relatively strong control measures should be given to groups with medium degrees. Furthermore, the allocation of costs (or resources) should coincide with their contact patterns.


Introduction
Measles, spread by coughing and sneezing, is one of the leading infectious diseases that causes death of children around the world. Vaccination is the key public health strategy to prevent it outbreaks. In 1986, China began to implement two doses of routine vaccination, at 8 months and 7 years old of age [1]. In 2010, a nationwide supplementary immunization activity (SIA) of measles was conducted in mainland China and subsequently the incidence of measles reached a low level of about 0.46/100,000 individuals in 2012 [2]. However, in 2013, the incidence of measles rose again in mainland China. Then the repeated outbreaks occurred with a pattern that was small-scale and clustered. Until now, the elimination goal that the incidence of measles is less than 1/100,000 population is still not achieved. Hence, how to design an effective strategy to eliminate the measles is still a challenge in mainland China.
Measles, like many other infectious diseases, exhibits great heterogeneity in terms of the number and the pattern of contacts. Much research indicated that heterogeneity may lead to highly different transmission dynamics of infectious diseases [3][4][5][6][7][8]. For example, the basic reproduction number, the threshold value to determine whether disease infection dies out or not, depends crucially on the variance of the distribution of contacts [9]. Accounting for the heterogeneity and non-random mixing, the meta-population basic reproduction numbers were about 70% greater [10]. Compared with the homogeneous massaction model, the extinction probability of the disease before a large epidemic occurs is much lower for the network model [11]. These results mean that the structure of the contact network is crucial in investigating and characterizing transmission of epidemics. So, many researchers studied the spread of epidemics such as foot-and-mouth disease, SARS, dengue, A/H1N1, and AIDS [12][13][14][15][16] on complex networks and analyzed the threshold of outbreaks and stability of equilibriums [17][18][19][20][21], which can seize the epidemic dynamics and further devise control measures. However, few studies have cautioned on the importance of the contact structures in understanding measles spread. Hence, it is necessary to investigate the network structure and how it relates to the propagation of measles. Moreover, for measles, waning of immunity may be of great risk to induce more infections, since vaccinated individuals with sufficiently low antibody levels are either at risk of disease infection [22] or have newborn babies with low immunity, who have potential risk to be infected before 8 months. Hence, how the waning of immunity does affect the measles infection also falls within the scope of this study.
Furthermore, an important topic of epidemiology, control and elimination of infectious diseases, have been extensively investigated for measles [23][24][25][26]. Particularly, in 1957, Bartlett [27] found that the localized extinction of measles was related to the community size and then it was observed by Keeling and Grenfell [28] in 1997 that this threshold may depend on the spatial structure and connectedness of the population. In 1999, Roberts and Tobias [29] found that the second MMR (measles-mumps-rubella) immunization at or around school entry may offer logistic advantages in New Zealand. In 2015, Verguet and Johri [30] proposed that without second routine immunization, regular SIAs at high coverage are able to control measles transmission, but the periodicity of SIA campaigns is determined by population demographics and existing MCV1 (measles-containing vaccine) coverage. Moreover, Hao and Glasser [31] calculated the gradient of the effective number with respect to age-and location-specific immunization rates and obtained young adults are the optimal target for SIAs. It should be noted that the majority of this literature for control of measles is based on the homogeneous mixed assumption. Hence, under heterogeneous scenario, how to design effective control measures for different population to eliminate measles is another aim of this study.
Our purpose of this study is to propose a novel network model with periodic transmission by dividing the vaccinated individuals into two classes: some vaccinated individuals with high antibody titer who are completely protected, others with low antibody titer who may be infected. We investigate the effective control strategy in this high-dimension and heterogenous system by applying optimal control method. We define the basic reproduction number and theoretically analyze the threshold dynamics of the system. Then we examine the optimal control problem and numerically investigate the effect of network structure on the outbreak of measles. In particular, we combine three kinds of feasible control measures (self-protection, enhancement of vaccination and treatment) and apply optimal control theory to answer the question of how to design, with the variation of network structure, the optimal control measures for each subpopulation with different activities and how much each group cost.

The network model
Each individual in a community can be regarded as a vertex in the network and each contact between two individuals is represented as an edge. The degree of a vertex means the number of contacts an individual has. In the following, the population is divided into n groups with sizes N k such that in group k, each individual has exactly k contacts for unit time. Assuming the population size is N , then the degree distribution of this network is defined as p(k) = N k N .

Model formulation
In this model, each vertex is empty or occupied by at most one individual. New recruitment with rate of constant b can only happen at the empty vertex and an occupied vertex becomes empty after the individual dies with rate of μ. Let S k , V 1k , V 2k , E k , E vk , I k , R k , 1 -A k be the densities of susceptible, vaccinated, exposed (infected from S k and V 2k ), infected, recovered and empty nodes with degree k. The vaccinated individuals (V 1k ) who are completely protected with high antibody titer progress to the class (V 2k ) with waning of immunity at rate of ω and then can be infected because of low antibody titer. We can write where Here p is the vaccination rate of the susceptible individuals, 1 σ is the incubation period of measles, γ is the recover rate of the infected individuals, the periodic function β(t) with period T is the baseline periodic transmission rate, η is the adjust factor which represents vaccinated individuals having reduced transmission rate (0 < η < 1). A flow diagram of the model is described in Fig. 1.  Based on the method provided in [32], we define the basic reproduction number of the system (1) as follows. Let y = (y 1 , y 2 , . . . , y 3M , y 3M+1 , . . . , y 7M ) represent the state of nodes in each class with y i ≥ 0 and divide the compartments into two types: infected classes

The threshold dynamics
. . , 3M and uninfected compartments (S k , V 1k , V 2k , R k ) relabeled by i = 3M + 1, . . . , 7M. Then we write the system (1) Denote , Specifically, Let Z(t, s), t ≥ s, denote the evolution operator of the following system: Then Z(t, s) satisfies where I is the identity matrix and the monodromy matrix -D (t) of the system (2) is equal to Z(t, 0), t ≥ 0. Assume that B T is the ordered Banach space composed of all T-periodic functions from R to R 3M with maximum norm · and then define a linear operator L : B T → B T as follows: Here, x ∈ B T denotes the initial distribution of infected individuals, G(s)x(s) represents the distribution of those newly infected by the infectious persons who were introduced at time s. For t ≥ s, Z(t, s)G(s)x(s) is the distribution of infected individuals produced at time s and remain in this compartment at time t. We can define the basic reproduction number of the system (1) by the spectral radius of L, i.e. R 0 = ρ(L).
For system (1), we can verify that the conditions (A1)-(A7) given in [32] hold and then derive the following theorem.

Theorem 1
The disease-free equilibrium E 0 is globally asymptotically stable if R 0 < 1, and unstable if R 0 > 1.
Proof Firstly, we illustrate that the disease-free equilibrium is locally stable. The Jacobian matrix of the system (1) at E 0 equals where with diagonal matricesÃ 1 ,Ã 2 ,B 1 ,B 2 ,B 3 , whose elements are -b, γ , p, -ω, -μ, separately. Clearly, all eigenvalues of J 2 have negative real parts. Then, combining with Lemma 1, that is, ρ( G-D (T)) < 1 for R 0 < 1, we derive that the disease-free equilibrium E 0 is locally stable if R 0 < 1. Next, we will illustrate E 0 is globally attractive.
For the system (1), it satisfies Then, using the standard comparison theorem [33], we can prove that, for any > 0, there Furthermore, we introduce an auxiliary system withˆ = m p(m|k)Î m and write it as and If R 0 < 1, corresponding to ρ( G-D (T)) < 1, we get ρ( G-D+ Q (T)) < 1 for enough small > 0 and consequently φ 1 = 1 T lnρ( G-D+ Q (T)) < 0. By Lemma 2.1 of [34], it can be illustrated that there is a positive T-periodic function c(t) = (c 11 (t), c 21 is a solution of the system (4). We can choose a small value α 1 > 0 and τ 2 > τ 1 to satisfyẐ(τ 2 Using the standard comparison theorem again, we get which demonstrates lim t→∞ E k (t) = 0, lim t→∞ E vk (t) = 0, lim t→∞ I k (t) = 0. Furthermore, considering the corresponding limit system Above all, we have completely proved the global attractivity of E 0 .

Theorem 2 If R 0 > 1, the system (1) is uniformly persistent and there exists at least a positive T-periodic solution.
Proof In the following, we first examine the uniform persistence of the system (1). This system admits a positive constantε such that, for all initial values with S 0 is the solution of the system (1) with initial value x 0 . Moreover, from the equations of system (1), we obtain lim t→∞ A k (t) = b b+μ , meaning that the solution of system (1) is uniformly and ultimately bounded. Hence, the semilow P is point dissipative and compact on R 7M + . Theorem 3.4.8 of [35] indicates that the map P has a global attractor, which attracts each , then by the equations of the system (1), we deduce for sufficiently small t > 0, v(t, x 0 ) / ∈ ∂Y 0 . This contradicts the initial assumption. Denote = ω(y 0 ), y 0 ∈ Y ∂ .
It is easy to see = {E 0 }, that is, for any initial value from Y ∂ , the solution of system (1) will remain in Y ∂ .
To show E 0 is a weak repeller for Y 0 , we only need to prove W s (E 0 ) ∩ Y 0 = ∅, where W s (E 0 ) is the stable manifold of E 0 . First, based on the continuity of solutions in terms of initial values, we see that ∀ > 0, there isδ > 0 such that, for any In the following, we demonstrate lim n→∞ sup d(P n (x 0 ), E 0 ) ≥δ.
Suppose not, there is some x 0 ∈ Y 0 such that lim n→∞ sup d(P n (x 0 ), E 0 ) <δ. Without loss of generality, we assume for any n > 0, Then the fourth and fifth equations of the system (1) satisfy Introducing the auxiliary system with˜ = m p(m|k)Ĩ m and for convenience, we rewrite (5) as ) and G(t), D(t), Q(t) have the same form as mentioned in Theorem 1. Following Lemma 2.1 [34], we can get a positive T-periodic functionc(t) = (c 11 (t),c 21 (t), , there exists a enough small > 0 such that ρ( G-D-Q (T)) > 1, which corresponds to φ 2 > 0. Choosingα 1 > 0, t 0 > 0 such thatZ(t 0 ) ≥α 1c (0), then the inequalityZ(t) ≥α 1c (tt 0 )e φ 2 (t-t 0 ) holds. Furthermore, the standard comparison theorem [33] demonstrates that the following inequality is established: This indicates lim m→∞ E k (t) = +∞, lim m→∞ E vk (t) = +∞, lim m→∞ I k (t) = +∞, which con- Above all, based on the Theorem 3.11 of [36], it is proved that the Poincaré map P is uniformly persistent with respect to (Y 0 , ∂Y 0 ). Furthermore, combining Theorem 1.3.6 [36], it can be proved that the Poincaré map has a fixed point , hence we can derive that the solution v(t, E * ) is a positive and T-periodic solution of the system (1). This completes the proof.

Numerical study
In this section, we initially analyze the effects of waning immunity on the transmission of measles in the random network and scale-free (SF) network. For the random network, it follows from Fig. 2(A) that increasing the waning rate (ω) of immunity can result in a significant increase in the basic reproduction number R 0 , indicating that waning of immunity brings about a difficulty in eliminating the measles. And it is worth noting that, given there is no waning of immunity, like in many existing models, one may underestimate the value of R 0 . Furthermore, Fig. 2(B) shows that an increasing waning rate leads to an increase in the density of infected nodes, especially the peak magnitude. This implies the waning of immunity causes many new infections and suggests us that prolonging the protection duration of vaccine is beneficial to control transmission of measles. Note that the plots for the SF network are similar to Fig. 2, and hence we omit them.
The basic reproduction number R 0 is the threshold to determine whether the measles dies out or not, hence it is essential to study the influence of the network structure determined by the degree distribution on R 0 . Here, we continue to discuss a scale-free (SF) network and a random network. In the SF network, the degree distribution follows the powerlaw distribution (that is, p(k) = mk -v , where m is uniquely determined by v) [37]. In the random network, the node degree follows the Poisson distribution (that is, p(k) = K k a e -Ka k! , the average degree K a = Nq, where N and q represent the network size and the probability of connection between two nodes) [38]. To study the impact of degree distribution on the basic reproduction number R 0 , we examine the variation in R 0 with parameter v (SF network) and q (random network) changing. Figure 3(A) shows that as the parameter v in power-law distribution increases, the proportion of nodes with large degree decreases, causing a reduction in R 0 . Hence, in the SF network, treating and isolating the infected individuals who have more connections to others can effectively control the spread of disease. Figure 3(B) indicates that in a random network, decreasing the probability q of random contacts between nodes can reduce the new infections.
To further explore the influence of network structure on the outbreak of disease, we compare the basic reproduction number R 0 between the SF network and the random network for a given average degree. Figure 3(C) indicates that there is a critical levelK a for the average degree (hereK a = 3.9) below which the basic reproduction number R 0 for the random network is greater than those for the SF network, above which the opposite result is observed. This means that the contact structure being like random network promotes the outbreak of measles when the average degree of nodes is small, while as the average degree increases, the proportion of nodes with large degree also increases and then measles can transmit faster in the SF network than in the random network.

Optimal control strategies
In this section, we consider some feasible interventions and then extend the system (1) to the following model with control functions u 1k (t), u 2k (t), u 3k (t): where W 1k = W 1 N k , W 2k = W 2 N k , W 3k = W 3 N k are the weights in terms of the control functions u 1k , u 2k , u 3k with positive constants W i , i = 1, 2, 3, k = 1, . . . , M. We can derive the optimal control strategies by searching for the optimal control function U * (t) subject to the system (6) such that . . , M} with positive upper boundŝ u ik , i = 1, 2, 3.

Optimal control and optimal solutions
Since the state system is a linear function for U and satisfies the Lipschitz property for the state variables, while the integrand of O is a convex function in terms of (u 11 (t), u 21 (t), u 31 (t), . . . , u 1M (t), u 2M (t), u 3M (t)), we can easily prove the existence of optimal controls by Corollary 4.1 of Fleming's reference [39] and derive the following result.

Theorem 3
There are the optimal control function U * (t) and the corresponding optimal solution C * (t) to minimize O(U(t)) over the˜ . In order to make the above statement correct, it is necessary to have continuous functions λ ik (t) such that with the transversality conditions, Furthermore, the optimal controls are u * 1k (t) = min max , 0 ,û 3k , k = 1, 2, . . . , M.

(8)
Proof Based on the Pontryagin's maximum principles [40], we can get the necessary conditions satisfied by the optimal controls. Specifically, we first define the Hamilton function H as The corresponding adjoint equations with transversality conditions are λ 1k = -∂H Differentiating H concerning u 1k , u 2k , u 3k on the set˜ , respectively, the Hamilton function H can reach the minimum with the following optimal conditions: Combining with the bounds 0 ≤ u 1k (t) ≤û 1k , 0 ≤ u 2k ≤û 2k , 0 ≤ u 2k ≤û 3k , k = 1, . . . , M, we derive the properties as Theorem 3.

Numerical simulation
In the following simulations, we mainly discuss the optimal control strategies and their corresponding effects on the incidence under two kinds of networks: random network and SF network. Here, we assume the community size N = 300 and the maximum degree M = 21.

Optimal control strategy on random network
To study the effect of optimal control strategies on the transmission of measles under different groups, we plot the time series of optimal control functions and the corresponding densities of infected nodes with degree k = 2, 11, 21, respectively. It is observed from Fig. 4(A)-(C) that the optimal control function u * 1k (t) begins to decline after a period of plateau with the highest level, and u * 2k (t) monotonically decreases while the optimal control function u * 3k (t) initially increases and then decreases which behaves like infection incidence. This demonstrates that we should initiate personal protection and enhance vaccination as soon as possible when disease takes off while the implementation of treatment measures should be in sync with transmission of measles. Moreover, the optimal control curves show that the nodes with high degrees require relatively long duration of self protection with maximum level and strong initial vaccination strength to minimize the objective function. Figure 4(D) indicates that with optimal controls, the outbreak can be significantly mitigated. Furthermore, the larger the degree k is, the higher the peak value of I k (t) is, which means the population with high degree is easily and seriously invaded by disease, so we should pay more attention to the individuals in these groups.
To further explore what interventions should be taken for different groups and their effectiveness on disease spread, we define the average control strength under the optimal controls and compare the cumulative densities of infected individuals I k ,Ī * k over a period under two cases: without/with optimal control, wherē It follows from Fig. 4 (A) that the average control strengthū ik increases with the degree k increasing. This indicates that the control measures related to the highly connected individuals should be strengthened to achieve optimal control. And further, for a given degree k, we haveū 1k >ū 2k >ū 3k , which means that the control strategies such as personal protection and social distance, associated with reducing transmission, should be significantly implemented for each group. Figure 5(B) demonstrates that the implementation of optimal controls yields a reduction of more than 50% in the cumulative density of infected individuals for each group.
To investigate the allocation of costs (or resources) in different groups under the optimal controls, we compute the fraction F i,k of costs occupied by each group for each strategy, as shown in Fig. 6(A), where  The distribution of each control costs among total groups in the random network (A) and SF network (B), which is illustrated by the fraction It displays that the cost distribution among different groups represents the shape of normal distribution and a substantially large fraction of the costs is occupied by the nodes with medium degrees, indicating that when the measles spreads among individuals whose contact structure is more like random network, more resources such as facial masks, vaccines, medicines should be allocated to the population with medium degrees. , we see that as the time delay increases, the duration that the optimal control u * 1,10 (t) being at the maximum level reduces while the initial intensity of optimal control u * 2,10 (t) also decreases. In contrast, the initial intensity of optimal control u * 3,10 (t) first increases and then declines with delaying the implementation of control strategies. Figure 7(D) displays that the delay in carrying out interventions leads to a significant increase in the densities of infected individuals. This indicates that the earlier the control measures is strengthened, the lower the outbreak size is.
Furthermore, we plot the average control strengthū 1k ,ū 2k ,ū 3k under four different starting times of control measures. Figure 8 shows that as the starting time is delayed, the average control strength of interventionsū 1k ,ū 2k declines, especially for the group with high degree k. This is because the relatively late implementation of control measures leads to the disease spreading out and then susceptible individuals quickly declining. In particular, the subpopulation with great connectivity is infected faster. Hence, the highest intensity of control measures associated with personal protection and enhancement of vaccination is implemented in the groups with medium degrees rather than in the groups with high degrees. This means that with time delay increasing, the control intensity for groups with medium degrees should be relatively stronger, compared to other classes, to achieve more effective control.

Optimal control strategy on SF network
To investigate the effect of different network structures on the control of disease, we compare the main results on scale-free network with those on random network. In contrast to the random network, we find from Fig. 6 (B) that, for the SF network, the costs among total groups follows the exponential distribution and a substantially large fraction of the costs is occupied by the nodes with low degrees, which suggests in SF network, more resources should be allocated to these corresponding groups. This indicates that the distribution of costs among various groups coincides with the distribution of contact pattern. In the SF network, most nodes have low degrees, which yields great resource requirement, while in the random network, most nodes are those with medium degrees and hence they occupy most resources.
Numerical studies exhibit that the optimal control functions in the SF network behave like those in the random network, and particularly, the nodes with high degree require high intensity of interventions. The specific plots are like Fig. 4 and we omit them. Similarly, with the starting time of control measures being delayed, the optimal control intensity for total groups decreases, especially for the group with large k, while the corresponding densities of infected individuals increase. Again, the plots for the SF network are almost the same as Fig. 7 and Fig. 8, hence we omit them.

Discussion
In recent years, the repeated outbreaks of domestic measles were of small-scale and clustered patterns, which presented obvious heterogeneity and attracted the public concern to design effective control measures against measles. To investigate the influence of heterogeneity on the measles transmission dynamics, we extend the homogeneous measles model to the heterogeneous network model. We further consider the waning of immunity in the proposed model and then divide the vaccinated individuals into two classes (fully protected with high antibody titer, probably infected with low antibody titer). We explore what control strategies should be implemented for individuals with different connectivities by using the optimal control method on the network model (high dimension). This network modeling approach provides a natural description to the transmission of measles among heterogeneous communities and we can obtain the detailed control scheme for individuals with different activities.
We initially analyze the threshold dynamics for the network model by defining the basic reproduction number R 0 , which is the spectral radius of the next generation operator. We prove that the disease-free equilibrium E 0 is globally asymptomatically stable for R 0 < 1. While, if R 0 > 1, the disease-free equilibrium E 0 is unstable and the system is uniformly persistent. The basic reproduction number R 0 can act as a threshold to determine whether the disease dies out or not. Then we numerically study the effect of waning immunity and find the larger the waning rate is, the greater the basic reproduction number R 0 and the density of infected individuals are, which demonstrates that the long-term protection of vaccine is effective in controlling transmission of measles. We investigate the influence of network structure on R 0 for two kinds of networks: the SF network (with degree distribution p(k) = mk -v ) and the random network (with degree distribution p(k) = K k a e -Ka k! , K a = Nq). Our results demonstrate that R 0 varies with the degree distribution for each network. The larger the parameter v or the smaller the parameter q of degree distribution is, the smaller R 0 is. This means we can effectively mitigate the outbreak of the disease by reducing the proportion of the highly connected individuals given contact structure in the case of the SF network, or restricting the contacts between nodes in the random network. Moreover, the comparison of R 0 between SF network and random network indicates that there exists a critical average degreeK a below which the basic reproduction number R 0 for the random network is larger than that, for the SF network, above which increasing the average degree induces a quick increase in R 0 (and thence new infections) for the SF network.
Furthermore, we extend the above network model by including three types of interventions (personal protection, enhancement of vaccination and treatment) to design the optimal control measures for achieving a low level of infections at a minimal cost. First, we prove the existence of the optimal solutions and characterize the optimal control functions using Pontryagin's maximal principle. The numerical results indicate that, no matter what the network is, the outbreak peaks early with relatively large magnitude for the subpopulation with high activity, hence these individuals require a high intensity of interventions.
With the optimal controls, the densities of infected nodes decrease by more than 50% for each group. Given delays of initiation of control measures, we see that late implementation of the interventions induces the relatively high density of infected individuals. This suggests that the control measures should be carried out as soon as possible. With the timing of implementation of control measures varying, the objective subpopulation implemented to the highest control strength also changes from individuals with high degrees to individuals with medium degrees, which is due to the delay of initiating controls results in the rapid spread of disease and further a rapid reduction of susceptible individuals, especially for those with great activities.
Our findings also demonstrate that, for the random network and the SF network, the costs occupied by different groups follow their own contact distribution. In particular, the distribution of costs is the normal distribution in the random network, while it follows the exponential distribution in the SF network. The results suggest us when the disease takes off in a community whose contact structure is more like random network, more resources should be allocated to the subpopulation with medium degrees due to most nodes are located in these groups. Conversely, in the SF network, most individuals have low degrees, so they should be paid more attention to and require more resources.
In summary, we developed a heterogeneous model with periodic transmission rate to investigate the effect of contact network and waning immunity on transmission of measles, which is an advantage compared to most existing homogeneous measles models. Furthermore, the effective control strategies for subpopulation with different activities are given by investigating the optimal control problem for the heterogeneous model. The results improve our understanding for the spread of measles and provide us with a detailed control scheme for different subpopulations, which cannot be obtained from the homogeneous models. However, it is challenging to obtain the realistic contact network due to the lack of reliable contact data, hence our simulations remain more qualitative. Moreover, it is essential to consider the effect of both age structure and contact heterogeneity on measles transmission, and we leave this for future work.