#### Abstract

This paper relies on the concept of next generation matrix defined ad hoc for a new proposed extended SEIR model referred to as -model to study its stability. The model includes successive stages of infectious subpopulations, each one acting at the exposed subpopulation of the next infectious stage in a cascade global disposal where each infectious population acts as the exposed subpopulation of the next infectious stage. The model also has internal delays which characterize the time intervals of the coupling of the susceptible dynamics with the infectious populations of the various cascade infectious stages. Since the susceptible subpopulation is common, and then unique, to all the infectious stages, its coupled dynamic action on each of those stages is modeled with an increasing delay as the infectious stage index increases from 1 to . The physical interpretation of the model is that the dynamics of the disease exhibits different stages in which the infectivity and the mortality rates vary as the individual numbers go through the process of recovery, each stage with a characteristic average time.

#### 1. Introduction

There is a very relevant interest in the literature concerning different aspects of dynamics of populations and related biological modeling issues including their positivity, stability, controllability, and observability and appropriate design of control rules for such models. The interest is twofold; namely, (a) on the one hand they have an undoubted mathematical interest because of the rich nonlinear dynamics they can exhibit which makes its analysis nontrivial in most cases, (b) on the other hand they have relevant interests in the real world concerning aspects such as health, resource exploitation, or rationalization of the labor management force towards economical issues. Biological and population dynamics models are usually defined by a variety of mathematical tools. Among them we can find differential and difference equations, eventually including internal delayed dynamics (i.e., delays in the state) or external delays (i.e., delays in the forcing action, if any) [1–8]. The delays can be, in general, modeled as point delays or as distributed delays and as constant or as time-varying delays. The background literature on the various involved subjects is exhaustive. In that context, important interest has been devoted to models of interaction of species versus their habitat, such as, for instance, the various Beverton-Holt models and some related generalizations, with their intrinsic problems of positivity, equilibrium analysis, stability, oscillatory solutions, and their control; see, for instance, [9, 10] and references therein. Important attention is paid to the study of their inverse models and their control which are equivalent to the initial direct ones while being much more tractable mathematically [9, 10]. Control rules through the online design of the habitat carrying capacity have been dealt with and applied in aquaculture exploitation [9]. Logistic equations, predator-prey models, and related oscillatory regimes have been also studied in the background literature (see, e.g., [11–13]). On the other hand, epidemic models of various kinds ranging from very elementary to sophisticated have received and still receive important attention. The literature is exhaustive (see, e.g., [1–8, 14–35]). Basic models are named via acronyms as SIR (susceptible-infectious-removed by immunity), SEIR (susceptible-infected-infectious-removed by immunity), or SVEIR (susceptible-vaccinated-infected-infectious-removed-by-immunity) models according to the involved subpopulations with coupled dynamics. An alternative popular name for the infected population is “the exposed” (those without external symptoms but already incubating the disease). An alternative popular name for the removed-by-immunity population is “the immune." There are also a wide variety of extended models available built with combinations and extensions of the above ones. The main properties dealt with are the calculation and related stability analysis of the disease-free and endemic equilibria, the infection permanence, which leads to the impossibility of reaching the disease-free equilibrium, and the positivity of the trajectories whose analysis combined with the boundedness of the total population leads to the internal stability [2, 3, 18, 20, 22, 24] of the whole model implying the boundedness of all the subpopulations [2–5, 18, 20, 22, 24–28, 30, 31, 35]. The control action on an epidemic model is performed through vaccination rules which can be of different types, for instance, constant or based either on linear or nonlinear feedback of some measurable or known subpopulations, mainly, the infectious population. Vaccination laws can also be of an impulsive nature (in practice, acting with large efforts along very short periods of time) or of a combined regular/impulsive nature [2, 5–7, 26–31]. In particular, the technique of adaptive sampling to design the relevant vaccination time instants is combined with the design of vaccination rules in [30] so as to increase the vaccination performance towards the disease eradication. It has to be pointed out that epidemic models are essentially uncontrollable so that it is not possible to drive simultaneously via a vaccination rule all the subpopulations to prefixed arbitrary values in finite prefixed time intervals since the total population is a constrain for all time for the sum of all the subpopulations [2, 7, 26]. However, most of the models are output controllable, or at least output stabilizable, with the outputs being defined “ad hoc” by either the infected and/or infectious subpopulations or the susceptible plus the immune subpopulations, those defined to be the output of the dynamic system, if necessary. The relevant idea in the epidemic models is that there is a transfer in-between subpopulations along the infection process (being, therefore, uncontrollable). The agent transmitting the infection is a quadratic dynamics of the susceptible and infected population with some either constant or functional factors (the incidence rate) which depends on a parameter or functional factor (the coefficient transmission rate) which depends on each infectious disease and the population nature [3, 16]. There is a very basic parameter to analyze, called the basic reproduction number, which has the following two interpretations.(a)If it ranges from zero to one, the disease-free equilibrium is asymptotically stable and the infection asymptotically vanishes without requiring any external action on the system. In this case, the Jacobian of the matrix of the dynamics of the epidemic model defining the linearized state-trajectory has all its eigenvalues in the stable region. If it exceeds one, the disease-free equilibrium is unstable so that the state trajectories can converge asymptotically to the endemic equilibrium or exhibit an oscillatory behavior.(b)The infection cannot progress if the reproduction number is less than one since the minimum available initial infected population to propagate the disease is below its critical value for spreading.

However, the explicit expression of the reproduction number is neither direct nor easy to obtain in sophisticated epidemic models. Therefore, more advanced techniques, based on the so-called next generation matrix, have been proposed to define such a reproduction number; see, for instance, [8, 33–35]. More recent studies on epidemic models are described in [36, 37]. The importance of demographics is taken into account in [36], while the stability of the equilibrium points of a delayed disease is described in [37]. This paper addresses the concept of the next generation matrix for a new proposed extended continuous time SEIR model referred to as model with delays where there are -successive stages of infectious subpopulations, each one acting at the exposed subpopulation of the next infectious stage in a cascade global disposal. Since the susceptible subpopulation is common to all infectious stages, its impact on each of those stages is modeled with an increasing delay as the infectious stage index increases from 1 to . A second related model is used without the above delay. Both models are also reformulated with discretized versions. The physical interpretation of the models is that the incubating period of any disease is not identical for all the susceptible individuals so that they can become infected at different times. So, each of the infectious stages is an infectious class which includes a certain number of individuals which become infected at times centered about a reference average time instant which is considered common for the whole group. The stability analysis is performed by defining a basic reproduction number from an ad hoc next generation matrix for this model. The basic properties which are proved are as follows. All the trajectories remain bounded since the differential system is nonnegative, as it is the alternative discrete model, in the sense that all the subpopulations are nonnegative for all time, and the total population is uniformly bounded for all time. If the defined basic reproduction number is less than one, then the disease-free equilibrium point is globally asymptotically stable as it is proved through the definition and use of a Lyapunov function. In this case, the endemic equilibrium point is locally unstable and, furthermore, it is unfeasible (i.e., nonreachable by any trajectory) since it is not compatible with the system positivity. If such a basic reproduction number is equal to one, then the endemic equilibrium does not exist as such since it is identical to the disease-free equilibrium and then globally stable since the system is positive and the total population is uniformly bounded for all the time so that none of the subpopulation can be unbounded. If the basic reproduction number exceeds one then the disease-free equilibrium point is locally unstable.

#### 2. The Model

This section contains the description and main properties of the introduced models. Previous models of multistaged infectious diseases have been described in [38] including (a) the so-called epidemic models with multistaged infectious period where the susceptible population influences the first infectious stage, that one, the second infectious stage, and so on; (b) the epidemic models with several types of infective subpopulations which are influenced in a parallel disposal and all of them have transitions to a unique removed population; (c) the epidemic among a number of homogeneous groups (each susceptible group generates its own SIR model). It can be considered, in a general context, that the structure proposed in this paper lies in the first of the above classes of the multistaged models.

##### 2.1. Model with No Delays

The proposed model is described as a succession of infectious stages of the disease, each one with a characteristic infectivity rate . The susceptible subpopulation goes through each infected subpopulation until they reach the recovered subpopulation, immune to the disease for a certain time until they become susceptible again. The susceptible population is born at a rate while the death rate of each infectious subpopulation is , with being the natural death rate of the population and the mortality rate associated with the disease at each infectious stage. The transition rate from the th infectious subpopulation to the next one as well as the transition rate from the last infectious subpopulation to the recovered one is denoted by , while the transition from the recovered to the susceptible is . All the parameters are assumed to be positive so as to possess full physical meaning. The equations for the dynamics of the subpopulations are coupled as follows: where and are the susceptible and recovered subpopulations and are the various infectious subpopulations corresponding to the different infection stages. For notation simplicity, the vector is defined as

Now the positivity, existence of equilibrium points, and stability are studied for this model.

###### 2.1.1. Positivity and Boundedness

Proposition 1. *All the subpopulations remain nonnegative for all time for any given nonnegative initial conditions, and all .*

*Proof. *Let be so that and . Then, from 1Now, let be so that and . Then, from 1For the next infectious subpopulations the same operation is made: let be so that and . Then, from 1Finally, let so that and , . Then, from 1
Since it has been proven that none of the subpopulations will have a negative derivative at any time instant when they vanish and the rest of the subpopulations are nonnegative, then all the subpopulations must have nonnegative values for all time.

Not only does the epidemic system possess nonnegative solutions under nonnegative initial conditions, but also the total population and all the subpopulations are bounded for all time with a uniform finite upper-bound as it is discussed in the subsequent two results.

Proposition 2. *The total population is finitely upper-bounded for any given nonnegative initial conditions.*

*Proof. *If we sum up 1, we obtain
Since Proposition 1 guarantees the nonnegativity of subpopulations and all the parameters are assumed to be positive, then as and for all time , for all and
The solution to the differential inequality 8 is given by
for all . Thus, the total population is finitely upper-bounded and the proposition is proved.

Corollary 3. *All the subpopulations are finitely upper-bounded for all time for any given nonnegative initial conditions, obeying the constraints , where is the total population at the disease-free equilibrium.*

###### 2.1.2. Equilibrium Points

The equilibrium points are obtained by zeroing the left-hand side of 1, resulting in This system possesses in general two solutions. The first one is the disease-free equilibrium point (DFE) in which the infectious and recovered subpopulations vanish so that the susceptible becomes the total population : and the endemic equilibrium point (END), described as is . Note that if , then , so the endemic equilibrium is not feasible so that the only reachable equilibrium point is the disease-free equilibrium point. In the next section, it is seen that this situation corresponds to the reproductive number to be less than one. Once we have the expressions of both equilibrium points, we are prepared to discuss their equilibrium in the next section.

###### 2.1.3. Stability Analysis: Next Generation Matrix

The local stability of the DFE point is proved in this section after introducing the reproductive number , which is the average number of new cases that produce an infected individual during the average duration of the disease. In order to find this number easily, a next generation matrix with small domain is constructed as follows. The Jacobi matrix defined around the disease-free equilibrium is composed of four different submatrices: , , , and :where would be the transmission matrix, representing the appearance of new infections, while represent the remaining transitions of the infected subpopulation unrelated as the deaths, or the transition from recovered to susceptible. Then, the inverse of that matrix would represent the average time in which a subpopulation stays in an infected state. The next generation matrix will be then defined as , and the reproductive number corresponds to the spectral radius of . Thus, the spectral radius is given by Thus, the stability of the disease-free equilibrium point is directly related to this spectral radius. In this way, it will be locally asymptotically stable when , while will imply instability of such a DFE point [35]. Studies of epidemic models with a parallel disposal of several submodels of susceptible-infectious-recovered subpopulations have proven that, for , the endemic point is locally asymptotically stable [39]. Remember that it has been proven in previous section, 12, that if , the only reachable point is the disease-free equilibrium, which is also locally asymptotically stable. If , both equilibrium points coincide (see 12).

##### 2.2. Discrete Model with No Delays

This section considers the discrete-time counterpart of the continuous-time model 1. In the discrete framework notation the subpopulations are written as for any time . The discrete model will be based on the continuous one from 1, with a step size . Thus, it is defined as the change , with as in the continuous subpopulations from 2. Therefore, it is obtained for that

Proposition 4. *The discrete system described in 17 is nonnegative for any nonnegative initial conditions if the step size is small enough to satisfy
**
with , , and , where and is the total initial population (i.e., the total population at the initial time).*

*Proof. *Given a set of nonnegative subpopulations at the th-sample, from the second equation of 17, it can be deduced that
Then, if . The same method can be applied for guaranteeing the nonnegativity of . From the third and fourth equations at 17, it is deduced that
Then, if .

In order to guarantee the nonnegativity of the susceptible subpopulation, first a maximum value is set for . From 17
then , and it can be said that
So from the first equation at 17 it is known that
Then, given that , , it is defined that so that the susceptible subpopulation remains nonnegative if
Thus, for a step size , the system remains positive.

The equilibrium points are obtained by making in 17:

The solutions of the above system of equations are the same as in the continuous-time case and given by 11 and 12 for the disease-free and endemic equilibrium points, respectively. A next generation matrix approach for discrete models is made so that the Jacobian matrix evaluated at the disease-free equilibrium is as follows: The Jacobian takes the formThe Jacobian can be separated into four submatrices as follows: with the fertility submatrix , related to the new infections and represented as and the transition submatrix , with the same dimensions, related to the transition of the individuals between the different subpopulations resulting in

A simple calculation with the Jacobian matrix results in an equivalent decomposition of the whole linearization about the equilibrium system into two complete subsystems, one describing the infection progress while the other one expresses the disease-free subpopulation dynamics as follows: The term represents the fraction of the individuals from the th infected subpopulation that will survive and move to the th. Because of these demographic interpretations, from Perron-Frobenius theory on the maximum modulus [40], the parameters are set so that the column sum of and the maximum moduli of its eigenvalues are less than one; that is, , so as to exclude the case of an immortal population. Then, it can be established for the vector of initial infectious individuals that . Then, the next generation matrix can be defined as the sum of the infections ever produced by the infected individuals at for all the time they remain infectious, which would be at , at , and at ad infinitum: where represents the distribution of all infections accumulated during the lifespan of the infectious population: Then the basic reproduction number is defined as the spectral radius of the matrix . The spectral radius of the is defined as . It is proven in [40] that either or or , so the stability of the disease-free equilibrium is determined by the value of . In our model this reproduction number is obtained as

Thus, the disease-free equilibrium state will be locally stable when the reproduction number is less than 1 and unstable otherwise. From known previous calculations of the stability of similar models [35, 40], this decomposition technique for stability analysis can be compared to previous ones such as the study of the stability made in [35] for the SEIR model, which is a specific version of a continuous model in which , , and . From the definition in 33 it is obtained that , which agrees with the study of the eigenvalues of the Jacobian in a continuous-time model [26].

##### 2.3. Model with Delays

###### 2.3.1. Construction of the Model

A new model is proposed based on the previous one in which the transition rates between subpopulations are substituted by delays, implying that each individual must stay at each stage of the infection during a certain period of time (latent period) before recovery. Such time periods are defined as for the infectious stages of the disease and for the recovered state. The model also preserves the infectivity mechanism of the previous one, in which all the infectious subpopulations affect the infection of the first stage at different rates. For notational abbreviation in the subsequent exposition it is now defined as the function . Thus, the value of the infected subpopulation at a certain time is defined by the ratio of the people that became infected at time that has not passed to the infectious subpopulation yet. In order to calculate such a total amount, consider the probability with being the probability distribution of the transition plus the probability to stay alive after time , which would be equal to . An integration is made in order to obtain the value of the subpopulation at the instant resulting in Given the above integral function of the first stage of the infectious subpopulations, its derivative through time is easily obtained as The probability distribution of transition is defined as , meaning that the probability of transition of the subpopulation at the stage 1 of the disease to the next one is set to zero before it has passed a latent time interval and is equal to 1 after that time, so that 36 transforms to As the rate of death in the subpopulation corresponds to , it is inferred that the term corresponds to the rate of transition from to . Then, the value of the people that arrives at at would be , and at certain time would be written as . For any subpopulation , it can be written that

For the recovered subpopulation the same technique is used, so that the dynamic of over time can be written as with , , and , . The dynamic equations of the subpopulations are defined then as where, as in the previous model, is the constant birth rate, and the fully immune recovered subpopulation eventually becomes susceptible again after the period of time .

###### 2.3.2. Positivity and Boundedness

Proposition 5. *The model equation 40 is nonnegative for any initial nonnegative conditions. Thus,
**, which implies ; .*

*Proof. *Since , it is then deduced from 41 that . Assume that
then, from 40 at the time ,
therefore, the subpopulation will not be the first subpopulation to be negative.

From 38 assume that such that for the first time and the value as small as desired such that ; . Then it can be said that
Since and knowing that is continuous and derivable for all , from the mean value theorem, it can be deduced that
which is in direct contradiction with the original assumption of ; . Then as is not negative before , the term will remain nonnegative. Consider . Therefore, neither the infectious nor the susceptible subpopulation will be the first one to have a negative value. The same method can be applied to demonstrate that the recovered subpopulation will not be the first to have a negative value either. Then, holds for all time for any given nonnegative initial conditions.

Since the lower bound of the subpopulations is established as 0, an upper-bound for all the subpopulations is set in the following proposition.

Proposition 6. *, , , for any given initial conditions .*

*Proof. *The dynamic of is obtained from 40: . As , and given the previous boundary of nonnegativity from Proposition 4, it can be seen that, for any , . Thus, . Then, for each subpopulation
it is established that ; .

###### 2.3.3. Equilibrium Points

The existence of equilibrium points in our model is now studied. The dynamics of the subpopulation described in 40 is zero at the equilibrium points, so that being . Two equilibrium points are obtained and result in with and . Note that and , so . Given any nonnegative initial conditions, the endemic equilibrium point would not be reachable if , as and would be <0.

Take into account that if , , the endemic, and the disease-free equilibrium are the same, so there is only one equilibrium point in the system.

###### 2.3.4. Stability of the Disease-Free Equilibrium (DFE) Point

The local stability of the DFE point is proved in this section after introducing the reproductive number , the average number of new cases that produce an infected individual during the average duration of the disease. In order to find this number easily, a next generation is constructed as follows.

First, the vector of the subpopulations is reorganized as in 2, so that the Jacobian matrix around the DFE point can be properly written as where is the transmission matrix, representing the appearance of new infections , while represent the remaining transitions of the subpopulation which are unrelated directly to the infection, such as the deaths, or the transition from recovered subpopulation to susceptible subpopulation.

Then, the next generation matrix is given by . The submatrices composing the Jacobian areso the transmission and transition matrices are
, . The element is the number of new cases generated in the stage of the disease by one infected case that has just arrived to the stage of the disease . Then, the dominant eigenvalue of the next generation matrix would correspond to the reproduction number . However, it is possible in this case to construct a small domain next generation matrix [35] with a lower dimension than , from which it will be much easier to obtain this dominant eigenvalue. Since (a) is a matrix whose rows are linearly independent vectors spanning the rows of and (b) is a matrix whose columns are linearly independent vectors spanning the columns of **,** then , and the small domain next generation matrix will be defined as . In our case and , so
since , the dominant eigenvalue of would be . Note that implies that if , then , which is also the condition of reachability of the endemic equilibrium, as stated in the previous section.

Proposition 7. *The DFE is globally asymptotically stable if .*

*Proof. *A candidate for a Lyapunov function is defined as , having the auxiliary functions , , been defined as
now, is positive definite , since , , , and . To verify the positivity of , consider the definition of the infectious subpopulation from 38 and the probability transition function . It is then established that
This integral function can be reconfigured as
Given this definition of , the parameter is defined as
Since and , it is deduced that , . Then and . Thus the nonnegativity of the function is proven. The function is equal to zero for , . Then the derivative over time is calculated as
. For , so for and otherwise from 40. Thus the proposition is proved.

###### 2.3.5. Stability of the Endemic Equilibrium Point (END) Model with

For simple models, that is, models with only one or two infectious subpopulations, a study of the stability of their endemic points is also made. The model is linearized around the endemic point and the eigenvalues of the Jacobi matrix are obtained. The function is defined as with the Jacobi matrix from 40, defined at the endemic equilibrium point from 49:

The solutions of are respectively, with and . Since and , and are defined as negative. The third solution will be defined as negative if , which is satisfied for or, as is shown in Section 2.3.4, when .

###### 2.3.6. Stability of the Endemic Equilibrium Point (END) Model with

For (i.e., there are two infectious subpopulations) a function will be defined as in 59, this time with a Jacobi matrix from the susceptible, the recovered, and the two stages of the infectious subpopulations: with being . The first eigenvalue is trivially obtained as . The other eigenvalues will not be directly obtained, as it is only needed to demonstrate that their real part is negative in order to prove the local asymptotic stability of the model. The Routh-Hurwitz criterion [41] says that in order to have all the solutions of on the left half plane, the coefficients of must satisfy the following conditions: (1)(2). Both conditions can be satisfied once that the limits of the parameter are established. Given and then so that, for , it is deduced that the coefficients are positive: And the lower limit of the second condition is defined as so all the real parts of the eigenvalues are negative and the equilibrium endemic point is locally asymptotically stable as long as .

##### 2.4. Discrete Model with Delay

The dynamic equation 40 is approximated to , with being the vector of the subpopulations to obtain a discrete-time counterpart of the continuous-time model. The notation for any subpopulation at the time interval is established then as . Thus, the analogous discrete equations to the model from 17 are obtained: and , , with and being and the rounded half down value for , , respectively. Now, at the equilibrium point, the delays are not relevant, while the dynamic equation at the equilibrium can be rewritten as

The reproductive number will be obtained again with a next generation matrix method as in the previous section, so the value would be which is similar to the reproductive number from the continuous model, with the only difference of the constants instead of the previous . However, given their similar origins and values, the properties exhibited in both sets of constants will be the same, in the sense that such as . Given the same model with the same characteristics, the values of the reproduction number will be very similar too.

#### 3. Simulation

A set of Matlab simulations are made based on the models described in the previous section, in order to contrast the predictions on their respective endemic and disease-free equilibrium states. A large enough number of infectious subpopulations are chosen, , so the stability of the models will be tested regarding their reproductive numbers. A common initial condition is set to all the susceptible subpopulation in the disease-free equilibrium plus a 0.1 of that value of infected subpopulation at the first stage of the disease. For exposition and calculations convenience, the mortality rate and the birth rate will be both taken as equal to the inverse of the average lifespan of a human . The mortality rate of the three stages of the infection will be , , , and the average time which an individual spends in each stage will be 19, 29, and 61 days, respectively. The infectious ratios will be altered in order to obtain reproductive numbers below and above 1, but the ratios between them will be constant, as . When simulating the discrete models, an appropriate time step size is chosen in order to guarantee the nonnegativity of model with no delays. As for the discrete model with delays, given that nonnegativity of the subpopulations in the continuous model has been proved in Proposition 5, an algorithm is established during the simulation so that the discrete dynamics approximates to the continuous one in any critical positivity points. A loop in which the following point is evaluated is introduced, acting over the time step size. While any of the subpopulations of the next point present a negativity for any , then the step size is reduced , with being the th loop iteration, , and . Another loop is set after this resetting the time step size while the next point is nonnegative or the time step size is equal to or below the original value.

##### 3.1. Model with No Delays

Given the parameters previously commented, the transition ratios between the different subpopulations are obtained as the inverse of their average times during each stage. Thus, , , , and the final transition from recovered to susceptible again will be set as . Given a value of the transition rates the reproduction number will be set to 0.5 and 1.5, respectively. For the continuous model, two graphics representing the evolution of each subpopulation for both situations are presented in Figures 1 and 2 while, for the discrete model, a step time is set to 0.1 years and the simulations are run with the same parameters. The results of these simulations are presented in Figures 3 and 4.

##### 3.2. Model with Delays

The common parameters are set as in the previous cases, but this time the delayed times will be equal to their respective average times years, years, and years. The delay of the transition between the recovered and the susceptible subpopulations will be changed from 610/365 to 80/365 years in order to let the simulation achieve the simulation at a reasonable time. The values of the transmission rates will maintain their proportionality, but the range of the possible reproductive number is not so wide as in the previous case. Their values are changed to now be and , respectively. The evolution of the subpopulations is presented in Figures 5 and 6.

In the discrete model, the same time step as before is set to years. The evolution of the subpopulations for each reproductive number can be seen in Figures 7 and 8, respectively.

Observe that the subpopulations rapidly tend to the disease-free equilibrium in Figure 8. While it is seen that the discrete-time and continuous-time models reach the same final states when the conditions are equal, the irregularity of the dynamics in the discrete models over the continuous has also been noticed. This fact reveals that the discretization procedure along with the time step must be selected carefully prior to simulating the discrete-time model.

#### 4. Conclusions

A model of a disease spreading with multiple infectious stages has been studied in depth. The next generation techniques have been proven to be quite useful when operating with models with complex interactions as in the variations of the model presented in this work. The reproduction numbers obtained in each case have also proven to predict correctly the final stable state of the subpopulations when they are properly simulated. The simulations have also shown that the different discretizations of a continuous model can cause significant discrepancies in the resulting dynamics. An appropriate integration step should be considered taking this into account as well as other factors like the computation time available, or, if the model is being used in a controlled system, the available time for data acquisition.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This work was partially supported by the Spanish Ministry of Economy and Competitiveness through Grant no. DPI2012-30651, the Scholarship support BES-2010-035160 for Raul Nistal, the Basque Government, through Grant no. IE378-10, and by the University of the Basque Country through Grant nos. UFI 11/07, PES13/21, and GIU14/07. Raul Nistal is grateful to the doctorate program “Molecular biology and Biomedicine” of UPV/EHU for its support for him as a Ph.D. student. Finally, the authors are also grateful to the reviewers for their useful suggestions.