Mathematical modeling of the transmission dynamics of an infectious disease with vaccination and optimal control analysis

Introduction

Compartmental models that divide a population into classes such as susceptible, exposed, infectious and recovered individuals have been the basic tool for describing how an infectious disease spreads through a population since the earliest work on the subject, and remain the starting point for evaluating the likely effect of an intervention such as vaccination before it is deployed. A model of this kind is useful in two related ways: it identifies the disease's basic reproduction number, the threshold quantity that separates conditions under which an introduced infection dies out from conditions under which it spreads, and it can be embedded in an optimal control framework, in which a public-health intervention is represented as a time-dependent variable to be chosen so as to minimize a stated objective, typically some combination of disease burden and the cost of the intervention itself. Adding an imperfectly effective vaccine to a susceptible-exposed-infectious-recovered structure, so that vaccinated individuals can still be infected, though generally at a reduced rate, is a standard extension, and computing the resulting reproduction number by the nextgeneration matrix method is now the standard approach for this class of model, with a general theorem establishing that the disease-free state is locally stable exactly when this reproduction number is below one. Treating the vaccination rate as a timedependent control rather than as a fixed constant, to be optimised via Pontryagin's Maximum Principle is likewise well established, and the resulting two-point boundary value problem-the original state equations forward in time coupled to a costate system that must be solved backward in time from a terminal condition-is routinely solved numerically by the forward-backward sweep method. This paper works through this standard pipeline end to end for a susceptible-vaccinated-exposed-infectious-recovered model, with two aims. First, every analytical step-the reproduction number, the stability threshold, and the costate equations arising from the Pontryagin Hamiltonian-is derived and then checked against an independent computation, rather than stated without verification, since it is easy for a sign or a factor to be lost when a costate system is derived by hand. Second, the resulting optimally controlled epidemic is compared numerically against the uncontrolled case for a disease whose reproduction number is high enough that, left alone, it would produce a substantial outbreak, to quantify how much of that outbreak a well-timed vaccination campaign can prevent.

Model Formulation

The population is divided into five classes: susceptible (S), vaccinated (V), exposed (E, infected but not yet infectious), infectious (I), and recovered (R). New individuals enter the susceptible class at a constant rate Λ (births or immigration) and leave any class through natural death at rate μ; infection adds an additional death rate d to the infectious class. Susceptible individuals are vaccinated at rate φ (a constant baseline rate in the reproduction-number and stability analysis, and a time-dependent control u(t) in the optimal control problem below) and move to the vaccinated class; vaccination is imperfect, reducing but not necessarily eliminating the risk of infection, captured by an efficacy parameter ε∈[0,1], so that a fully effective vaccine (ε=1) prevents infection entirely while ε=0 corresponds to no protection. New infections occur at rate β per contact with an infectious individual, scaled by the fraction of the population that is infectious, and move susceptible or vaccinated individuals into the exposed class at rates β and (1-ε)β respectively. Exposed individuals progress to the infectious class at rate σ (the reciprocal of the mean incubation period), and infectious individuals recover at rate γ. With N=S+V+E+I+R denoting the total population, the model is

Basic Reproduction Number and Stability Analysis

At the disease-free equilibrium (E=I = 0), setting dS/dt=dV/dt= 0 gives S*=Λ/(μ+φ), V*=Λφ/[μ(μ+φ)], and N*=S*+V*=Λ/μ. Following the next-generation matrix method, with the exposed and infectious classes ordered (E,I), the new-infection matrix F and the transition matrix 𝒱, evaluated at the disease-free equilibrium, are

and the basic reproduction number is the dominant eigenvalue of F𝒱⁻¹,

This expression was checked in two ways. Setting φ=0 (no vaccination) reduces it exactly to the classical susceptibleexposed-infectious-recovered reproduction number σβ/[(μ+σ)(μ+γ+d)], confirming that the vaccination-related terms vanish correctly when there is no vaccination. Setting ε=1 (a perfectly effective vaccine) and letting φ→∞ (universal, immediate vaccination) drives R0 to zero, confirming that complete, perfect vaccination eliminates transmission in the model, as it must. The general threshold theorem for compartmental disease models states that the disease-free equilibrium is locally asymptotically stable if R0 < 1 and unstable if R0 > 1; for the present model, the full fivedimensional Jacobian at the disease-free equilibrium has two eigenvalues, -μ-φ and -μ (the latter with multiplicity two), that are negative for any admissible parameters, together with a pair of eigenvalues governing the exposed-infectious subsystem whose sign matches the sign of 1-R0 by the general theorem. This was confirmed directly by evaluating the full Jacobian numerically at two parameter sets: one giving R0 = 0.142, for which all five eigenvalues were negative, and one giving R0 = 2.564, for which one eigenvalue was positive, in both cases matching the stability status implied by R0.

Optimal Control Formulation

The baseline vaccination rate φ is now replaced by a timedependent control u(t), bounded between 0 and a maximum feasible rate umax, and chosen to minimise the objective functional J(u) = ∫₀ᵀ [ A1 E(t) + A2 I(t) + (B/2) u(t)²] dt, over a fixed control horizon [0,T], where A1 and A2 weight the cost of exposed and infectious individuals respectively and B weights the cost of the vaccination effort itself; the quadratic control cost is the standard choice in this literature and reflects diminishing returns to a faster vaccination campaign. For tractability, the total population N in the transmission terms is treated as a constant over the control horizon during the derivation of the costate equations, a standard simplifying assumption in this literature; the full state equations, with N left to evolve as S+V+E+I+R, are used without this simplification in the numerical solution reported below. Forming the Hamiltonian with costate (adjoint) variables λS, λV, λE, λI, λR attached to the five state equations and applying Pontryagin class does not feed back into any other equation and does not appear in the objective functional, its costate equation is a linear homogeneous equation with a zero terminal condition, whose unique solution is λR(t)≡0 for all t; this is a genuine simplification of the problem rather than an approximation. Minimising the Hamiltonian with respect to u and applying the bounds gives the characterisation of the optimal control, u*(t) = min{ umax, max{ 0, S(t)[λS(t)-λV(t)]/B } }.

All of the derivatives above were obtained and independently re-derived using symbolic computation to guard against sign or algebra errors before being used in the numerical scheme described next.

Numerical Method

The coupled two-point boundary value problem-state equations with initial conditions at t=0, costate equations with terminal conditions at t=T-was solved by the forwardbackward sweep method. Starting from an initial guess for u(t) (zero, in the results below), the state equations are integrated forward in time with a fourth-order Runge-Kutta scheme; the costate equations are then integrated backward in time from the terminal condition, using the just-computed state trajectory; the control is updated at every time point from the optimality characterisation above, using a convex combination of the new and previous control values for numerical stability; and the process repeats until the largest change in the control across the whole time grid falls below a fixed tolerance (10⁻⁶ in the results below, reached after 68 iterations on a grid of 1601 points over 160 days).

Results

The following baseline parameters were used for the optimal control simulation, chosen to represent a moderately transmissible disease with a multi-day incubation and infectious period: Λ=2 per day, μ=1/(70×365) per day, β=0.5 per day, ε=0.85, σ=1/5 per day, γ=1/7 per day, d = 0.01 per day, with objective weights A1 = 1, A2 = 3, B = 60000, and a maximum vaccination rate umax = 0.15 per day over a 160day horizon, starting from S(0)=9990, V(0)=0, E(0)=5, I(0)=5, R(0)=0. In the absence of vaccination, these parameters give a basic reproduction number R0 = 3.270, indicating a disease capable of substantial spread if unchecked.

Table 1: Comparison of key outcomes with no vaccination versus the optimal vaccination control

No vaccination

Optimal control Peak number infected 1870.5 (day 55.0) 13.9 (day 22. Figure 1 shows the exposed and infectious trajectories on a logarithmic scale, since the two scenarios differ by more than two orders of magnitude at their peaks. Without vaccination, the infectious population grows for about eight weeks before peaking at around 1870 individuals and then declining as susceptibles are depleted; with the optimal control in place, the infectious population never exceeds about 14 individuals, and both the exposed and infectious curves stay close to their initial values throughout. Figure 2 shows the optimal vaccination rate itself. It starts just below the imposed upper bound umax = 0.15 per day, without ever reaching that bound exactly, and decreases smoothly and monotonically over the following weeks as the risk of a large outbreak recedes, reaching a negligible level by around day 100. The optimal strategy is therefore front-loaded: nearly all of the vaccination effort recommended by the model is concentrated in the first few weeks, before exponential growth in the number of infectious individuals has had time to take hold, rather than being spread evenly across the horizon or deployed reactively after infections have already risen.

Discussion

The scale of the difference between the controlled and uncontrolled outcomes in Table 1 is a direct consequence of the exponential nature of early epidemic growth: with R0 = 3.270, every infectious individual who is prevented early in the outbreak also prevents the several further generations of infection that individual would otherwise have caused, so a comparatively modest, front-loaded vaccination effort (an average rate well below the 0.15 per day upper bound, applied mainly in the first three to four weeks) is far more effective than the same total effort would be if it were spread evenly across the full horizon or delayed until after the infectious population had already begun to grow quickly. This is consistent with the general finding, common across this literature, that the timing of an intervention relative to the exponential growth phase of an outbreak matters at least as much as its total magnitude. The stability analysis and the optimal control analysis are complementary rather than independent results: the reproduction number R0=3.270 computed for the uncontrolled system explains why the epidemic in Table 1 grows as sharply as it does in the first place, and it is exactly this threshold quantity that the vaccination control must, in effect, drive below one, locally in time, for the outbreak to be suppressed rather than merely slowed. The two internal checks reported in Sections 3 and 4-the correct reduction of R0 to the classical formula and to zero under perfect universal vaccination, the numerically confirmed match between the sign of the Jacobian eigenvalues and the R0 threshold, and the analytical reduction of the recovered-class costate to zero-were included specifically to guard against the kind of algebraic slip that a hand derivation of a five-compartment optimal control problem can easily introduce, and all three checks were consistent with the underlying theory before the model was used to generate the results in Table 1 and Figures 1 and2.

Limitations

The model assumes homogeneous mixing within a single population, exponentially distributed sojourn times in each compartment, and a vaccine whose effect is captured by a single efficacy parameter rather than separate parameters for susceptibility and onward transmission. The optimal control problem uses a constant population size when deriving the costate equations, while the numerical state simulation lets the population vary through births, deaths, and disease mortality; for diseases with substantial mortality or over very long horizons this simplification would need to be revisited. Only vaccination is treated as a control; other interventions such as treatment, isolation, or contact reduction, which are often modelled alongside vaccination in this literature, are not included. Finally, the quadratic cost of control and the specific weights A1, A2 and B are illustrative choices; a real application would calibrate these to the costs and constraints of an actual vaccination programme, and the qualitative frontloading result reported here, while consistent with the underlying exponential-growth mechanism, would need to be re-checked against the specific parameters of the disease and programme in question.

Conclusion

This paper formulated a susceptible-vaccinated-exposedinfectious-recovered model with imperfect vaccination, derived its basic reproduction number and disease-free stability threshold, and set up and solved an optimal control problem for the vaccination rate using Pontryagin's Maximum Principle and the forward-backward sweep method, checking each analytical step against an independent computation along the way. For a disease with a basic reproduction number of about 3.3 in the absence of vaccination, the optimal, frontloaded vaccination strategy reduced the peak infectious population by more than two orders of magnitude and cumulative incidence by a factor of more than seventy relative to no vaccination, reflecting the outsized value of intervening before exponential growth takes hold rather than after. The same verified pipeline-next-generation matrix, stability check, Pontryagin derivation, and forward-backward sweepextends directly to models with additional compartments or additional controls, such as treatment or isolation, without requiring a different analytical approach.

Acknowledgement

AcknowledgementsThe author(s) thank the Department of Public Health Mathematics, Gulf Medical Research University, Doha, Qatar for providing the computational resources used in this work. All symbolic and numerical computations were carried out using open-source scientific computing libraries (SymPy, NumPy and Matplotlib).

References

1. Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London Series A. 1927;115(772):700721.
2. van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences. 2002;180:29-48.
3. Diekmann O, Heesterbeek JAP, Roberts MG. The construction of next-generation matrices for compartmental epidemic models. Journal of the Royal Society Interface. 2010;7(47):873-885.
4. Pontryagin LS, Boltyanskii VG, Gamkrelidze RV, Mishchenko EF. The Mathematical Theory of Optimal Processes. New York: John Wiley & Sons; 1962.
5. Lenhart S, Workman JT. Optimal Control Applied to Biological Models. Boca Raton: CRC Press; 2007.