Critical behavior in a stochastic model of vector mediated epidemics

The extreme vulnerability of humans to new and old pathogens is constantly highlighted by unbound outbreaks of epidemics. This vulnerability is both direct, producing illness in humans (dengue, malaria), and also indirect, affecting its supplies (bird and swine flu, Pierce disease, and olive quick decline syndrome). In most cases, the pathogens responsible for an illness spread through vectors. In general, disease evolution may be an uncontrollable propagation or a transient outbreak with limited diffusion. This depends on the physiological parameters of hosts and vectors (susceptibility to the illness, virulence, chronicity of the disease, lifetime of the vectors, etc.). In this perspective and with these motivations, we analyzed a stochastic lattice model able to capture the critical behavior of such epidemics over a limited time horizon and with a finite amount of resources. The model exhibits a critical line of transition that separates spreading and non-spreading phases. The critical line is studied with new analytical methods and direct simulations. Critical exponents are found to be the same as those of dynamical percolation.

VME by coupled SIR-SIS models. In order to describe the spreading of VME, we adopt the recently introduced coupled SIR-SIS model 10 , which we simply refer to as the VME model. The SIR and SIS dynamics interact in the spreading mechanism, like in real life. This model captures most of the relevant mechanisms of epidemic spreading, although in a schematic way. The VME model involves two different players, hosts (H) and vectors (V). Only a cross mechanism of infection is available (catalytic infection). Hosts follow SIR evolution, while vectors are governed by SIS dynamics. This is inspired by the idea of two species with very different life expectancies. Hosts and vectors are placed on a regular checkerboard (see 20 for epidemics in complex networks). Self-interactions are absent, and the evolution of one subspecies is always mediated by another species in neighborhood. Susceptible players of each species are infected with a specific virulence rate. The model is schematically illustrated in Fig. 1. As we briefly mentioned in the introduction, most studies point out the global-warming-induced increase in vector populations and in pathogen virulence as one of the main reasons for recent epidemic outbreaks and interest (see, for instance [21][22][23][24]. The SIS component present in our formulation implies instead a one-to-one reproduction of vector units. Thus, our modeling would apply to a local, closed population within a limited time horizon and with limited vital resources available (such as water, or any other reproduction/incubation media).
The different rates of infection for hosts and vectors are denoted by v H and v V respectively, which generally differ. The effect of host vaccination is implemented by reducing the rate v H . Vectors are assumed to die more rapidly than hosts. After being infected, they do not recover but die and are reintroduced in the game as susceptible players, with exit rate e V . The effect of prophylaxis is to increase this rate and reduce the number of infected vectors. Finally, the infected hosts may be put out of the game (recover or die) with an assigned exit rate, e H . For a small exit rate, hosts tend to remain in the infected state, and their illness becomes chronic. In contrast, a large exit rate reduces the number of individuals that they can infect. Note that the rate of exit may also be associated with a prophylactic practice, like the eradication of infected trees in VMEs like OQDS 25 . Tuning these ratios up or down leads to winning of the spreading or extinguishing phases. A list of the aforementioned rates and of some other useful symbols is reported in Table 1.

Infection paradigms.
To reduce the number of parameters, we adopt the condition v H − e H = v V − e V , which represents the equality of the effective rates of spreading for each species 10 . There is no deep epidemiological reason behind this condition that is a mere simplifying assumption. In principle, our analysis could be extended to the more involved 3-parameter case. Anticipating our findings, we shall see that each point of the critical line will be equivalent from the point of view of critical behavior being characterized by the same exponents governing its size dependence -although with varying critical amplitudes. Switching on more parameters, as it would be in a realistic model, replaces the critical line by a critical surface where the same features are expected. In our setup, the independent variables reduce to the ratios e H /v V , v H /e V , and the epidemic spreading can be explored in the plane of p = 1/(1 + e H /v V ) and h = 1/(1 + v H /e V ). For each specific ratio of the host exit and vector virulence rate, it is possible to find a corresponding ratio of the vaccination and prophylaxis rates sufficient to counteract the spreading. Conversely, the non-spreading region NS is reached for each value of h, making the ratio e H /v V larger than its critical value. In a realistic context, this view would be more significant in VMEs with non-human hosts, such as OQDS or Pierce Disease, for which vaccination is still not allowed and improving the host exit rate is an admissible prophylactic strategy. In other words, in order to confine the epidemics, it is enough to make the exit rate (eradication) just a little bit larger than (1/p c − 1)v V .

Results
The critical line. The continuous set of critical points (p c , h c ) constitutes the critical line of the model (see where P(Λ ; t) is the probability that the full-lattice VME model is in the state Λ at time t, and w(Λ → Λ ′ ) are the rates of the configuration change Λ → Λ ′ . Assuming invariance under translations, the Master Equation can be turned into an infinite coupled system of evolution equations for the joint probabilities P(s 1 , s 2 , … ) for finding a cluster of neighbouring sites in the states s i . The stochastic model is thus substantially beyond any simple description that does not take into account fluctuations on the microscopic scale. In 10 , an approximate critical line has been determined by a stability analysis of the Master Equation. This required a closure hypothesis to deal with a finite number of evolution equations. The mean-field (MF) approximation corresponds to the uncorrelated choice P(s i , s j ) = P(s i ) P(s j ). The pair mean field (PMF) line is obtained by replacing the three-site probabilities with the combination = P s s s ( , , ) i j k P s s P s s P s In the MF case, the evolution of the system is truncated to a set of three coupled differential equations, while in the PMF case gives a system of five equations. For the analytical determination of the critical line, see Methods A. In panel (a) of Fig. 2 it is represented by the yellow (MF) and red (PMF) curves. Neither approximation is accurate in reproducing the data points (green circles) obtained by Monte Carlo simulation 10 . We have further improved this approximation scheme by introducing two arbitrary numerical coef- The factors τ i,k are regarded as phenomenological parameters that encode the effect of the higher joint probabilities. In principle, they may be determined by comparison with Monte Carlo simulation data 26 . The best improved critical line is obtained by the choice τ = .
1 438 and is shown by the blue curve in Fig. 2. Remarkably, all data points are accurately fit by our simple two-parameter modification.
Panels (b-d) of Fig. 2 give a qualitative idea of the general spatio-temporal characteristics of the VME epidemics in different regions of the (p, h) plane. This is illustrated by snapshots of typical realizations of the epidemics on a 1000 × 1000 lattice. Susceptible hosts and vectors are represented by light gray and white pixels, while infected hosts and vectors are represented by red and orange pixels, respectively. In panel (b), we consider a point near the critical line but in the non-spreading phase. We show the cluster of recovered sites (RSC) at the time where the epidemics stop. The cluster has a fractal structure and does not reach the lattice boundary. In panel (c), we are precisely on the transition line. The RSC is yet again fractal and shown at an arbitrary intermediate time. At later times, the periodic conditions push the cluster to wrap around the boundaries until the epidemic stops. Finally, in panel (d), we show a realization associated with a point well within the spreading phase. Again, we show the RSC at an intermediate time. The cluster is now compact with an approximately circular shape. The geometrical structure of the infected sites changes as we move off the critical line. At criticality in panel (c), there are small disconnected clusters of infected hosts and vectors at the boundary of the growing cluster of recovered hosts. Instead, deep in the spreading phase in panel (d), these clusters merge in a connected ring along the boundary. A quantitative separation of the spreading and non-spreading phases requires a detailed finite size scaling analysis, to be given later, see equation (3). For more details about the temporal evolution, see Fig. 3. Additional information related to the reproductive ratio Finite Size Scaling data analysis and critical exponents. The evolution of the VME model is investigated by using a Monte Carlo simulation on a checkerboard square lattice with L × L sites, assuming periodic boundary conditions. We start with one infected host at the centre of the lattice and let the infection evolve until there are no more infected sites. This is an inactive (adsorbing) state. For each realization of this stochastic process, we measure various quantities of interest. The order parameter P is defined to be one if the infection reaches the lattice boundary and zero otherwise. The number N H,R is the total number of recovered hosts when the infection stops (which always happens in finite volume). These static quantities are averaged over all realizations, leading to 〈 P〉 and 〈 N H,R 〉 which are functions of L, p, h. Following 18 , we also introduce the important combination 2 . This is analogous to Binder's cumulants 27 and may be shown to be scale invariant. During the simulation, we also record the number of infected hosts N H,I when a configuration reaches the boundary, along with the elapsed time T at that moment. The averages 〈 N H,I 〉 and 〈 T〉 are computed by restricting to realizations with a spreading infection. We explore the critical line by fixing h and varying p. For each value of h, we expect to find a critical p c (h) where the spreading-non-spreading transition occurs. From Finite Size Scaling FSS theory (e.g., 28 ), we expect the following laws: where ε = p − p c (h), β, γ, and ν ⊥ are universal critical exponents, and the coefficients k P , ′ k P , k H,R , and ′ k H,R are non-universal metric factors that are expected to be dependent on the spectator coupling h 29 . F P and F H,R are universal finite-size scaling functions in the sense that their dependence on L, p, h is fully encoded in the indicated variables. A similar relation holds for the cumulant B, but since B is scale invariant, it has no L-dependent pre-factor and its FSS relation is: where the value of B at ε = 0 is also universal. The meaning of equations (3) and (4) is rather simple. For an infinite system, the space correlation length ξ diverges like ξ ε ν − ⊥ . In finite size L, the free energy depends on the coupling through the ratio L/ξ up to a global power of L providing the right dimension of the observable quantity under consideration. In general, equations (3) and (4) are empirically valid in a wide region around the transition point ε = 0, but scaling corrections are typically present, although subdominant for large L. To determine p c (h) and the exponents, we exploit equation (3) because, for ε = 0, it predicts a simple linear behavior between the logarithm of the scaling quantity and log L. An alternative method to locate p c (h) is based on equation (4). It implies that the curves plotting B as a function of p at different values of L intersect precisely at the point ε = 0.
As mentioned, we analyze the critical line by evaluating the (constant) exponents in equations (3) and (4) at different values of h with a suitable choice of the (technically non-universal, i.e. model dependent) metric constants k and k′ .
We have run the VME model on lattices with sizes up to 1024 and a typical number of 10 7 realizations. The analysis has been repeated at the four values h = 0.1, 0, 2, 0.4, 0.6, i.e., a growing rate of vaccination for an assigned rate of prophylaxis, or equivalently, a growing rate of prophylaxis for an assigned rate of vaccination. In Fig. 4, we show the linear behavior of the quantities log〈 P〉 and log〈 N H,R 〉 plotted as functions of log L for the four choices of h. Our best fit results for the critical points and exponents are reported in Table 2. The universal exponents are the same along the critical line. Their values enable us to conclude that, concerning the exponents, the present VME model belongs to the same class of the SIR model 18 . Precisely, we find values of β/ν ⊥ and γ/ν ⊥ that are in good agreement with the known values β = 5 36 , γ = for Dynamical Percolation (DynP) 15 . The same agreement is also found for the cumulant B, which was precisely measured to be B = 1.0167(1) 18 . A second test is shown in Fig. 5, where we show the collapses of the data for the four values of h, according to equations (3) and (4). For any value of h the data collapse in a smooth curve, confirming the scaling hypothesis. Moreover, with a suitable choice of the metric factors, the four curves overlay precisely, giving strong support to the universality of the three functions F P , F H,R , and F B . Comments about the dynamical critical exponents can be found in Methods B.

Discussion
We have improved the determination of the critical line exhibited by the VME model proposed in 10 and identified the character of the transition along the whole line, showing that the model has the critical exponents as standard dynamical percolation. In our setup, the two rate ratios p and h play similar roles, and each point corresponds to a specific disease or specific prophylactic practices. Hence, the asymmetric shape of the critical line may suggest that to improve the epidemic confinement, improving h could be more efficient rather than reducing p. To give an interpretation, this means that prophylactic practices acting on the vectors (the rate e V ) may be more efficient than prophylactic practices acting on hosts, like eradication of trees in the case of the OQDS or Pierce Disease infection. Of course, to make concrete these remarks one should consider more realistic and deep modelizations that are not limited to the analysis of the critical dependence on finite resources, being just one facet of the problem. Then, it should become necessary to account for the actual costs of the different counteractions. Besides, in a real context, ethical problems are a major issue and must be taken into account. From this point of view, any suggestion for a realistic public health policy will definitely separate human infections and agricultural ones, since ethical issues and admissible interventions in the two cases are completely different. As we stressed, this is out of reach for our VME model. Nevertheless, from the point of view of the critical behavior, we have shown the irrelevance of the SIS sub-dynamics. The VME model is actually a SIR model in disguise. Similar conclusions are expected to hold for analogous models appearing in completely different applications, like in the context of internet viral diffusion 30 , opinion dynamics 31 , or traffic simulation 32 .  (p, h) plane. We show three snapshots taken during the temporal evolution of the VME model on a 1000 × 1000 lattice at the three points p = 0.76, 0.78, 0.9 (rows from top to bottom), and h = 0.2, as considered in Fig. 2. The first row is in the non-spreading phase. After the initial outbreak, the number of infected hosts initially grows and then reduces, leading to a final static state where the cluster of recovered hosts has not reached the lattice boundary. The second row is on the critical line. The infection invades the lattice, wraps around the periodic boundary, and eventually terminates. The final cluster of recovered hosts is fractal due to isles of different shapes containing untouched susceptible hosts surrounded by recovered ones. Finally, the third row shows the evolution deep in the spreading phase. Now, the evolution leaves behind a compact cluster of recovered hosts. In the final state, the lattice is almost filled with such sites.  According to the scaling forms in equations (3) and (4), the critical exponents β/ν ⊥ and γ/ν ⊥ may be extracted from the slope of the data in log-log scale and are clearly independent on the value of the coupling h. The error bars in the second plot are smaller than the dots.  probability P(I V ) for a vector to be infected decreases according to the exit rate for vectors e V , and increases through contacts with infected hosts. The gain term is proportional to the probability of finding the vector in the susceptible state, with a first neighbouring infected host. The interpretation of the other two lines in equation (5) is similar. If we include the next level in the hierarchy, i.e. the equations for the 2-sites probabilities, the system becomes H V V H 0 A study of the explicit numerical solutions of the system (6) shows that, as expected, all epidemics end in a finite time in both phases, reaching a stationary states without infected hosts or vectors. The temporal evolution of R 0 (t) is non trivial. In the non-spreading phase, R 0 (t) drops monotonically to a constant asymptotic value < 1, and its first derivative is negative and monotonic as well. On the other side, in the spreading phase, R 0 (t) shows an almost constant plateau, with R 0 (t) > 1, during the initial spreading of the epidemic, and eventually drops as soon as the stationary state is reached. The fact that all epidemics can be eradicated within a finite time horizon is strictly due to the fact that we simulate in finite size, i.e. with a local closed population. This would be an unsound assumption for the realistic all-round study of widespread and global endemic infections such as malaria, dengue and various others. In the present context, this is not a problem because we are only focusing on the study of the size dependence and need to extrapolate from finite size data. The study of the fluctuations of the quantities appearing in equation (10) in the lattice stochastic model is an interesting issue that could provide improved estimators for the spreading/non-spreading transition, but is beyond the scope of the present study.
Dynamical exponents. Apart from the static critical exponents, the DynP class is also characterized by a further dynamical exponent. This is defined at criticality in infinite size L → ∞ by the temporal evolution of the growing lattice of recovered hosts. This dynamical exponent appears in the scaling law 〈 T〉 ~ L z , where z = ν ‖ /ν ⊥ measures the ratio of correlation exponents in time and space [33][34][35] . A precise measure of z may be found in 36  is the fractal dimension of the recovered hosts cluster. Our analysis is focused on the static, i.e. spatial, critical properties and is less precise concerning dynamical quantities. From our data for 〈 T〉 , we obtain z = 1.13(1) and d F − z = 0.76 (9) for the scaling of 〈 N H,I 〉 , in fair agreement with the expected values 37 .