Deterministic Seirs Epidemic Model for Modeling Vital Dynamics, Vaccinations, and Temporary Immunity

In this paper, the author proposes a new SEIRS model that generalizes several classical deterministic epidemic models (e.g., SIR and SIS and SEIR and SEIRS) involving the relationships between the susceptible S, exposed E, infected I, and recovered R individuals for understanding the proliferation of infectious diseases. As a way to incorporate the most important features of the previous models under the assumption of homogeneous mixing (mass-action principle) of the individuals in the population N, the SEIRS model utilizes vital dynamics with unequal birth and death rates, vaccinations for newborns and non-newborns, and temporary immunity. In order to determine the equilibrium points, namely the disease-free and endemic equilibrium points, and study their local stability behaviors, the SEIRS model is rescaled with the total time-varying population and analyzed according to its epidemic condition R0 for two cases of no epidemic (R0 ≤ 1) and epidemic (R0 > 1) using the time-series and phase portraits of the susceptible s, exposed e, infected i, and recovered r individuals. Based on the experimental results using a set of arbitrarily-defined parameters for horizontal transmission of the infectious diseases, the proportional population of the SEIRS model consisted primarily of the recovered r (0.7–0.9) individuals and susceptible s (0.0–0.1) individuals (epidemic) and recovered r (0.9) individuals with only a small proportional population for the susceptible s (0.1) individuals (no epidemic). Overall, the initial conditions for the susceptible s, exposed e, infected i, and recovered r individuals reached the corresponding equilibrium point for local stability: no epidemic (DFE XDFE) and epidemic (EE XEE).


Introduction
Over the past decades, mathematical models have been developed and implemented to study the spread of infectious diseases since the early 20th century in the field of mathematical epidemiology [1][2][3][4]. The stochastic and deterministic epidemic models allow researchers to gain valuable insights into numerous infectious diseases and investigate strategies for combating them. For deterministic epidemic models, the populations of individuals are assigned to one of several different compartments, where each compartment represents a specific stage of the epidemic. The transition rates from one compartment to another compartment are mathematically expressed with derivatives. Based on the assorted compartments for the population and derivatives for the transition rates, the system of ordinary differential equations serves to describe the changes in population as a function of time.
From the seminal work in 1927, Kermack and McKendrick constructed a simple deterministic compartment model that still today acts as the fundamental model for developing and implementing even more complicated mathematical epidemic models [5]. In Kermack and McKendrick's SIR classical compartmental model, the population N was divided into three compartments: susceptible S compartment, in which all the individuals are susceptible if they have contact with a disease; infected I compartment, in which all the individuals are infected by the disease and infectious to spread the disease; and recovered R compartment, in which all the individuals are recovered from the infection. For the SIR model, Kermack and McKendrick made three assumptions: (1) the disease spreads in a closed environment (i.e., no births or deaths) with constant population N; (2) the number of susceptibles S who are infected by an infected I individual per unit of time is proportional to the total number of susceptibles with the proportional (transmission) coefficient β, where total number of newly infectives is βSI; and (3) the number of recovered R individuals from the infected I compartment per unit time is γI, where γ is the recovery rate coefficient, with the recovered individuals gaining permanent immunity from the infectious disease. While the SIR model was accurate for describing the spread of viral diseases (e.g., influenza, measles, and chickenpox), the SIR model was inappropriate for dealing with bacterial diseases (e.g., encephalitis and gonorrhea), where the recovered individuals gained no immunity and could be re-infected by the disease again at a future time period. As a spinoff of the SIR model, Kermack and McKendrick proposed the SIS model 5 years later to study the dynamics of bacterial diseases [6]. In both the classical SIR and SIS models, the models assume a negligible disease incubation period, where the susceptible S individuals could become infected and later recovery to acquire permanent or temporary immunity.
For more general models than SIR and SIS models, the SEIR and SEIRS models assume that the susceptible S individuals-after infection-first proceed through the latent period as exposed E individuals before becoming infected I individuals and then eventually recovered R individuals [7,8]. In the exposed E compartment, the individuals are infected by the disease but do not have the visible symptoms of the disease and cannot communicate the disease to susceptible S individuals. In such a latent period, the disease takes a certain time for the infection to multiple inside the body of the susceptible S individuals to reach the critical level to become infected I individuals. After the incubation period of the disease, the exposed E individuals soon become infected I individuals and then either acquire permanent immunity (SEIR) or temporary immunity (SEIRS). As with the SIR and SIS models, the SEIR and SEIRS models assume homogeneous mixing (mass-action principle) of the individuals in the population.
With the SIR and SIS models along with the SEIR and SEIRS models, the models all assumed that the disease spreads in a closed environment. For such models, the population N is always a constant value since the models do not incorporate any births or deaths. In order to develop and implement more realistic mathematical models to mimic reality, Anderson and May investigated the use of vital dynamics to vary the size of the population [9,10]. By assuming a birth rate b and death rate d, the SEIR and SEIRS models with vital dynamics now have a time-varying population N(t) that more appropriately models the spread of the disease. Fundamentally, the total population increases by birth at a rate b and decreases at a rate of d. From the different compartments of susceptible S compartment, exposed E compartment, infected I compartment, and recovered R compartment, the total population N(t) adheres to the conversation law as the sum of the populations of the different compartments that all vary as a function of time.
In the presence of infectious diseases, the ideal goal is to fully eradicate them through either preventive measures or establishment of a mass immunization program. As an extension to the SEIR and SEIRS models with vital dynamics, Anderson and May studied vaccinations applied to newborns (i.e., babies) and non-newborns (i.e., children and adults) [11,12]. For mass immunization programs, newborns or susceptible S individuals receive the vaccines and proceed directly to the recovered R compartment. By providing the proper vaccines to the public, the mass immunization program serves to reduce the basic reproduction value R 0 to less than unity, which causes the infectious disease to die out eventually. For R 0 greater than unity, the infectious disease does not die out eventually and actually causes the occurrence of an epidemic. Due to successful immunization programs occasionally creating health problems to individuals since vaccinations offer some associated risks, the infected I individuals may become too fearful of the risks and then not receive the necessary vaccines until they have spread the infectious disease to many other susceptible S individuals. Before reaching an epidemic, legislation is sometimes passed to enforce vaccinations.
As a way to merge the significant features of the foundational and more advanced SIR and SIRS models [13,14] along with the SEIR and SEIRS models [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34], the focus of the work is to develop and implement an extension of the SEIRS model. Fundamentally, the new SEIRS model is now a more advanced generalization of the previous models and incorporates vital dynamics with unequal birth and death rates, vaccinations for both newborns and non-newborns, and temporary immunity for describing the spread of infectious diseases. The SEIRS model with vital dynamics, vaccinations, and temporary immunity is rescaled using the total time-varying population and analyzed to determine its equilibrium points and corresponding local stabilities of the equilibrium points. In order to test the SEIRS model, numerical simulations are run involving a set of arbitrarily-defined parameters for horizontal transmission of the infectious disease in the new SEIRS model.
The remainder of this paper is organized into the following sections: epidemiological model (Section 2), local stability (Section 3), experimental methodology (Section 4), experimental results (Section 5), and conclusion (Section 6).

Epidemiological Model
In epidemiology, the new SEIRS model with vital dynamics (birth and death rates), vaccinations (newborns and non-newborns), and temporary immunity provides a mathematical description of infectious diseases and corresponding spread in biology. Figure 1 shows the block diagram of the SEIRS model with compartments (classes) consisting of susceptible S, exposed E, infected I, and recovered R individuals from the total population N.
occasionally creating health problems to individuals since vaccinations offer some associated risks, the infected I individuals may become too fearful of the risks and then not receive the necessary vaccines until they have spread the infectious disease to many other susceptible S individuals. Before reaching an epidemic, legislation is sometimes passed to enforce vaccinations.
As a way to merge the significant features of the foundational and more advanced SIR and SIRS models [13,14] along with the SEIR and SEIRS models [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34], the focus of the work is to develop and implement an extension of the SEIRS model. Fundamentally, the new SEIRS model is now a more advanced generalization of the previous models and incorporates vital dynamics with unequal birth and death rates, vaccinations for both newborns and non-newborns, and temporary immunity for describing the spread of infectious diseases. The SEIRS model with vital dynamics, vaccinations, and temporary immunity is rescaled using the total time-varying population and analyzed to determine its equilibrium points and corresponding local stabilities of the equilibrium points. In order to test the SEIRS model, numerical simulations are run involving a set of arbitrarily-defined parameters for horizontal transmission of the infectious disease in the new SEIRS model.
The remainder of this paper is organized into the following sections: epidemiological model (Section 2), local stability (Section 3), experimental methodology (Section 4), experimental results (Section 5), and conclusion (Section 6).

Epidemiological Model
In epidemiology, the new SEIRS model with vital dynamics (birth and death rates), vaccinations (newborns and non-newborns), and temporary immunity provides a mathematical description of infectious diseases and corresponding spread in biology. Figure 1 shows the block diagram of the SEIRS model with compartments (classes) consisting of susceptible S, exposed E, infected I, and recovered R individuals from the total population N.
The four classes of the model, S, E, I, R, are described in further detail in Table 1.  For the SEIRS model, the model allows for vital dynamics with unequal birth and death rates, vaccinations of both newborns and non-newborns, and temporary immunity from the infectious disease. With the population, the newborns are all susceptible without vaccination. In specific terms, the newborns are not born with any maternally-derived immunity. At a later time period, the newborns (i.e., babies) who were not initially vaccinated and non-newborns (i.e., children and adults) The four classes of the model, S, E, I, R, are described in further detail in Table 1. For the SEIRS model, the model allows for vital dynamics with unequal birth and death rates, vaccinations of both newborns and non-newborns, and temporary immunity from the infectious disease. With the population, the newborns are all susceptible without vaccination. In specific terms, the newborns are not born with any maternally-derived immunity. At a later time period, the newborns (i.e., babies) who were not initially vaccinated and non-newborns (i.e., children and adults) who were not around during the inception of the infectious disease as newborns both compromise the susceptible class and have the ability to receive vaccines. Based on the interaction of susceptible individuals and infected individuals, the infectious disease transmission causes the susceptible individuals to now leave the susceptible class and enter the exposed class. After a specific time period, the incubating infection eventually causes the exposed individuals to obtain the infectious disease and then spread it to susceptible individuals. Without any vaccines to temporarily cure the exposed individuals and infected individuals, the infected individuals will recover from the infectious disease and enter the recovered class. Unfortunately, the recovered individuals experience only temporary immunity from the infectious disease and can potentially transition back into the susceptible class. Table 2 summarizes the interpretations of the different positive parameters embedded in the SEIRS model for each of the four classes. Rate at which susceptible individuals become exposed by infected individuals and leave susceptible class and enter exposed class σ Transmission Rate (Exposed to Infected) days −1 Rate at which exposed individuals become infected by incubating infection and leave exposed class and enter infected class γ Recovery Rate (Infected to Recovered) days −1 Rate at which infected individuals leave infected class and enter recovered class From the positive parameters, the rates ρ (vaccination rate of non-newborns), α (transmission rate of recovered to susceptible), β (transmission rate of susceptible to infected), σ (transmission rate of exposed to infected), and γ (recovery rate of infected to recovered) are numerically interpreted in terms of their inverses. Table 3 explains the interpretations of the positive inverse parameters embedded in the SEIRS model. Based on the model parameters and inverse model parameters in Tables 2 and 3, the SEIRS model is transformed into mathematical system for analysis and evaluation. Mathematically, the SEIRS model is expressed as a system of ordinary differential equations given as: with population: or: where β SI N is the incidence rate at which the susceptible S become infected I by a disease. By substitution of the ordinary differential equation system in (1) into the relationship of (3), the population N is governed by the ordinary differential equation: Since (4) does not depend on any of the other variables in the system in (1), the population N is computed using separation of variables: or: to yield: As a result of solving (7) through exponentiation, the population N is given as: with time-varying population N(t).
Instead of solving the ordinary differential equation system in (1) with known population N from (8), the transformations: and: and: and: are applied to the system in (1), where s, e, i, and r denote the fractions of the number of individuals in classes S, E, I, and R with population N. Now, the transformed system (see Appendix A for more details) is formulated as: which is equivalent to the system in (1). By substitution of the transformations in (9)-(12), (2) is written as: or: With manipulation of (14) to produce: (16) is substituted into the transformed system in (13) to eliminate r and yield the simplified subsystem: or: with positive constants: and: Based on the solutions of the subsystem in (18), s, e, and i are utilized to solve for r in: is a positive constant. In order to transform the system in (13) back to the original system in (1), the solutions s, e, i, and r and N in (8) are inserted into the transformations in (9)-(12).

Local Stability
From the transformed subsystem in (18), the local stability is analyzed to determine the disease-free equilibrium (DFE): and endemic equilibrium (EE): Specifically, the equilibrium points are computed by setting ds dt = 0, de dt = 0, and di dt = 0 and solving for s, e, and i in (18) to compute the two equilibrium points. From the equilibrium points, the Jacobian matrix J is calculated from: as:

∂G/∂s ∂G/∂e ∂G/∂i ∂H/∂s ∂H/∂e ∂H/∂i ∂K/∂s ∂K/∂e ∂K/∂i
and evaluated at the equilibrium points to decide on the local stability, which is directly determined from the eigenvalues λ of: Based on the eigenvalues λ of (29), the linearized system will either be stable (all the eigenvalues of the Jacobian evaluated at the equilibrium point contain negative real parts) or unstable (at least one of the eigenvalues of the Jacobian evaluated at the equilibrium point has positive real part) for the transformed subsystem in (18).

Disease-Free Equilibrium
Through substitution of (25) into ds dt = 0 of the transformed subsystem in (18), the DFE X DFE in (25) is computed as: to generate: or: At the DFE X DFE = A B , 0, 0 , the Jacobian matrix J(X) in (28) is given as: with eigenvalues λ: J X DFE − λI = 0 (34) or: which is expanded as: By evaluating the determinants in (36), the eigenvalues λ are determined from the cubic polynomial: where the three eigenvalues λ are dependent on the parameters σ and β and constants A, B, C, and D.
Unfortunately, the eigenvalues λ are difficult to compute for the cubic polynomial in (37) without any specific values for the parameters σ and β and constants A, B, C, and D. In order to determine the parameter and constant independent local stability of the DFE X DFE in (25), the Routh-Hurwitz criteria is applied to the cubic polynomial in (37) with the coefficients: and: Based on the Routh-Hurwitz criteria for a cubic polynomial P(λ), the three conditions: and: must satisfied for the DFE X DFE in (25) to be locally stable. For the first condition (41): since the constants B > 0, C > 0, and D > 0 in (20)- (22). With the second condition (42): From the third condition (43): or: By cancellation of BCD and −σβA in (45) and (48), (43) is simplified as: or: After dividing by BC, (50) is simplified as: where: As a result of the Routh-Hurwitz criteria, all the eigenvalues λ in the cubic polynomial P(λ) in (37) have negative real parts to conclude that the DFE X DFE in (25) is locally stable with (46).

Endemic Equilibrium
From the transformed subsystem in (18), the EE X EE in (26) is computed by developing a relationship between i and e through di dt = 0 as: where: By examining de dt = 0 as: or: (55) is substituted into (57) to deliver the first coordinate of EE X EE as: With ds dt = 0 as: or: (58) is inserted into (60) as: After distributing and collecting terms in (61): which is further simplified with (54) as: or: By distributing and combining fractions, (64) is written as: or: to supply the second coordinate of EE X EE as: Through substitution of (67) into (54), the third coordinate of the EE X EE is given as: With (58), (67), and (68), the EE X EE in (26) is given by: which only makes physical sense if: since all the constants A, B, C, and D and parameters α, β, and σ in (69) are positive values. By manipulating the inequality in (70), the epidemic condition R 0 is given as: The epidemic condition R 0 in (71) is the basic reproduction value and is the most important quantity to consider for analyzing any epidemiological model. In particular, R 0 determines whether an epidemic occurs for infectious diseases since R 0 is the average number of secondary infections produced by one infected individual during the mean period of infection in a fully susceptible population. If R 0 ≤ 1, then, on average, the number of new infections produced by one infected individual over the mean course of the infectious disease is less than unity, which implies the infectious disease dies out eventually. Conversely, if R 0 > 1, then, on average, the number of new infections produced by one infected individual is greater than unity, which leads to the persistence of the infectious disease as an epidemic. At the EE X EE in (69), the Jacobian matrix J(X) is given by: which is expanded as: Through evaluation of the determinants in (75), the eigenvalues λ are computed from the cubic polynomial: or: where the three eigenvalues λ are dependent on the parameters σ and β, constants A, B, C, and D, and first and third coordinates of EE X EE , namely s* and i* in (69). In a similar manner to the eigenvalues λ for the cubic polynomial in (37), the eigenvalues λ for the cubic polynomial in (77) are even more difficult to compute without any specific values for the parameters σ and β and constants A, B, C, and D. The Routh-Hurwitz criteria with conditions (41)-(43) is again applied to now the cubic polynomial in (77) to determine the parameter and constant independent local stability of the EE X EE in (69) with the coefficients: and: From the first condition (41): since the constants A > 0, B > 0, C > 0, and D > 0 and parameters α >0, β > 0, σ > 0 in (20)- (22) and table. With the second condition (42): if (71). For the third condition (43): Through multiplication of (83) and (84) by (CD + ασ + αD) 2 , (43) is simplified as: or: As a consequence of the Routh-Hurwitz criteria, all the eigenvalues λ in the cubic polynomial P(λ) in (77) have negative real parts to conclude that the DFE X EE in (69) is locally stable with (71).

Experimental Methodology
The proposed SEIRS model with vital dynamics, vaccinations, temporary immunity developed in (1) and transformed to the subsystem of ordinary differential equations in (18) was evaluated in Matlab. Table 4 lists a typical set of numerical values for the model parameters and inverse model parameters for all the experiments. Based on the numerical values of the model parameters and inverse model parameters, the positive constants A, B, C, D and F ((19)- (22) and (24)) along with the epidemic condition R 0 ((71)) are computed for the SEIRS model in (18) and shown in Table 5. For the SEIRS model in (18) with rescaled variables s, e, i, and r ((9)-(12)), the total population n is assumed to have a value of unity ( (14)), where the population N of the original SEIRS model in (1) is assumed to have a value of 100 ((2)).

Experimental Results
Simulations were performed using the proposed and rescaled SEIRS model with vital dynamics, vaccinations, and temporary immunity in (18) along with the numerical values of the model parameters and inverse model parameters (Table 4) and constants (Table 5). In order to differentiate between the possibilities for the epidemic condition R 0 , the two cases of no epidemic (R 0 ≤ 1) and epidemic (R 0 > 1) were analyzed separately. For both the no epidemic and epidemic cases, the local stability of the DFE X DFE and EE X EE equilibrium points in (32) and (69) were evaluated using their corresponding eigenvalues λ i from their specific Jacobian matrix J X . The proportional population of the rescaled variables s, e, i, and r with total population n of unity was studied through the subsequent time-series with various initial conditions s(0), e(0), i(0) and r(0) over the course of a 90-day period. As a way to examine the relationships between the different rescaled variables, phase portraits were utilized to trace the solution of the system of ordinary differential equations for the rescaled SEIRS model in (18). With all of the numerical simulations, the time period is assumed to have units of days.

No Epidemic
Based on the numerical values of the model and inverse model parameters and constants with β = 1/4 (0.25 susceptible individuals who become exposed by infected individuals and leave the susceptible class and enter the exposed class per day) for the rescaled SEIRS model in (18), the epidemic condition R 0 was calculated as R 0 = 0.21865, which is less than unity and implies no epidemic for the infectious disease. Table 6 shows the DFE and EE equilibrium points X DFE and X EE and eigenvalues λ i of the Jacobian matrices J X DFE and J X EE along with their local stabilities. From the eigenvalues λ i of the DFE and EE equilibrium points X DFE and X EE , the DFE equilibrium point X DFE is locally stable (all eigenvalues with negative real parts) and the EE equilibrium point X EE is locally unstable (at least one eigenvalue with non-negative real part) for β = 1/4.  is locally stable (all eigenvalues with negative real parts) and the EE equilibrium point EE X is locally unstable (at least one eigenvalue with non-negative real part) for β = 1/4. Figure 2 illustrates the time-series of the proportional populations for the rescaled variables against time in days for a 90-day time period using various individual initial conditions.   With the assorted phase portraits, the initial conditions (s(0), e(0), i(0), r(0)) for the rescaled variables s, e, i, and r always reach the DFE equilibrium point DFE X since it is locally stable for the case of no epidemic (β = ¼). With the assorted phase portraits, the initial conditions (s(0), e(0), i(0), r(0)) for the rescaled variables s, e, i, and r always reach the DFE equilibrium point X DFE since it is locally stable for the case of no epidemic (β = 1 4 ).

Epidemic
From the numerical values of the model parameters and inverse model parameters and constants with β = 4 (4 susceptible individuals who become exposed by infected individuals and leave the susceptible class and enter the exposed class per day) for the rescaled SEIRS model in (18), the epidemic condition R 0 was calculated as R 0 = 3.49840, which is greater than unity and implies epidemic for the infectious disease. Table 7 displays the DFE and EE equilibrium points X DFE and X EE and eigenvalues λ i of the Jacobian matrices J X DFE and J X EE along with their local stabilities. In contrast to the no epidemic case (β = 1 4 ), the epidemic case (β = 4) causes the DFE and EE equilibrium points X DFE and X EE to have different locally stability. With the eigenvalues λ i of the DFE and EE equilibrium points X DFE and X EE , the DFE equilibrium point X DFE is locally unstable (at least one eigenvalue with non-negative real part) and the EE equilibrium point X EE is locally stable (all eigenvalues with negative real parts) for β = 4. Figure 4 exhibits the time-series of the proportional populations for the rescaled variables against time in days for a 90-day time period using various individual initial conditions.   Through the proportional population time-series, the exposed e and infected i individuals do not eventually decay to zero as with the no epidemic case (β = 1 4 ). In fact, exposed e and infected i individuals reach their steady-state of 0.0-0.2 in approximately 10-40 days. At the end of the studied time period of 90-days, the proportional population consists of recovered r (0.7-0.9) individuals and susceptible s (0.0-0.1) individuals. Figure 5 presents the two-dimensional phase portraits of the various combinations of the rescaled variables using the initial conditions of (s(0), e(0), i(0), r(0)) = (0.25, 0.25, 0.25, 0.25).   By the variety of phase portraits, the initial conditions (s(0), e(0), i(0), r(0)) for the rescaled variables s, e, i, and r always reach the EE equilibrium point EE X since it is locally stable for the case of epidemic (β = 4).

Conclusions
In this paper, the author developed and implemented a new SEIRS model that capitalized on the mutual benefits of the SIR and SIRS and SEIR and SEIRS models. Fundamentally, the focus was By the variety of phase portraits, the initial conditions (s(0), e(0), i(0), r(0)) for the rescaled variables s, e, i, and r always reach the EE equilibrium point X EE since it is locally stable for the case of epidemic (β = 4).

Conclusions
In this paper, the author developed and implemented a new SEIRS model that capitalized on the mutual benefits of the SIR and SIRS and SEIR and SEIRS models. Fundamentally, the focus was to generalize the previous models to incorporate vital dynamics with unequal birth and death rates, vaccinations for newborns and non-newborns, and temporary immunity for communicating the advancement of infectious diseases. From the experimental results of the proposed and rescaled SEIRS model, the local stability of the DFE X DFE and EE X EE equilibrium points were examined for the cases of no epidemic (R 0 ≤ 1) and epidemic (R 0 > 1) using the time-series and phase portraits of the susceptible s, exposed e, infected i, and recovered r individuals. Whereas the exposed e and infected i individuals eventually decayed to zero after approximately 40 days (no epidemic), the exposed e and infected i individuals reached their steady-state of 0.0-0.2 in approximately 10-40 days (epidemic). In the no epidemic case, the proportional population consisted of recovered r (0.9) individuals with only a small proportional population for the susceptible s (0.1) individuals. For the epidemic case, the proportional population consisted primarily of the recovered r (0.7-0.9) individuals and susceptible s (0.0-0.1) individuals. For both the no epidemic and epidemic cases, the initial conditions for the susceptible s, exposed e, infected i, and recovered r individuals reached the corresponding equilibrium point for local stability: no epidemic (DFE X DFE ) and epidemic (EE X EE ). For future work, the SEIRS model with vital dynamics, vaccinations, and temporary immunity could be modified to incorporate age structure, infection-age structure, and spatial structure along with treatment, isolation, quarantines, and vertical transmission to obtain even more realistic epidemic mathematical models.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A
Based on the transformations in (9)-(12), the SEIRS model in (1) with vital dynamics (birth and death rates), vaccinations (newborns and non-newborns), and temporary immunity is formulated mathematically as a system of ordinary differential equations as in (13). By substitution of (9)-(12) into dS dt , dE dt , dI dt , and dR dt of (1): and: and: and: dR dt S = sN ) or: After insertion of (4), (A5)-(A8) are rewritten as: Through cancellation of N and collecting terms, (A9)-(A12) are given as ds dt , de dt , di dt , and dr dt in the transformed ordinary differential equation system in (13).