# SIR Epidemic Model with Vaccination ## Overview This model simulates the spread of an infectious disease through a population using the SIR (Susceptible-Infected-Recovered) compartmental model with vaccination intervention. The model demonstrates how vaccination programs can reduce disease transmission and save lives. ## Mathematical Model The system is governed by the following coupled differential equations: ``` dS/dt = -βSI/N - νS (Susceptible) dI/dt = βSI/N - γI (Infected) dR/dt = γI + νS (Recovered/Removed) ``` Where: - `S(t)` = number of susceptible individuals at time t - `I(t)` = number of infected individuals at time t - `R(t)` = number of recovered/removed individuals at time t - `N = S + I + R` = total population (constant) ## Parameters ### Disease Parameters - **`β = 0.5`** - **Transmission rate** - Units: 1/time - Rate at which susceptible individuals become infected upon contact with infected individuals - Higher values indicate more contagious diseases - The term `βSI/N` represents new infections based on contact between S and I populations - **`γ = 0.1`** - **Recovery rate** - Units: 1/time - Rate at which infected individuals recover (or are removed from the population) - Reciprocal gives the average infectious period: `1/γ = 10` days - Higher values indicate faster recovery/removal ### Intervention Parameters - **`ν = 0.01`** - **Vaccination rate** - Units: 1/time - Rate at which susceptible individuals are vaccinated - The term `νS` represents people moving directly from S to R through vaccination - Higher values indicate more aggressive vaccination programs ### Population Parameters - **`N = 1000`** - **Total population** - Units: individuals - Constant total population size - Used for normalization in the transmission term `βSI/N` ## Key Epidemiological Metrics ### Basic Reproduction Number **`R₀ = β/γ = 0.5/0.1 = 5.0`** - Represents the average number of secondary infections caused by one infected individual in a fully susceptible population - `R₀ > 1`: Epidemic will spread - `R₀ < 1`: Epidemic will die out - `R₀ = 1`: Epidemic threshold ### Herd Immunity Threshold **`H = 1 - 1/R₀ = 1 - 1/5 = 0.8 = 80%`** - Fraction of population that must be immune to stop transmission - Can be achieved through natural infection or vaccination ## Initial Conditions ### Scenario 1: `iv1 = [990, 10, 0]` - 990 susceptible individuals - 10 initially infected individuals - 0 initially recovered individuals - Represents a small initial outbreak ### Scenario 2: `iv2 = [990, 20, 0]` - 990 susceptible individuals - 20 initially infected individuals (double the initial outbreak) - 0 initially recovered individuals - Used to demonstrate sensitivity to initial outbreak size ## Simulation Parameters - **Time span**: 0 to 200 days - **Step size**: `h = 0.1` days - **Numerical method**: 4th-order Runge-Kutta (vectorRK4) ## Model Variations ### With Vaccination ``` dS/dt = -βSI/N - νS dI/dt = βSI/N - γI dR/dt = γI + νS ``` ### Without Vaccination (Comparison) ``` dS/dt = -βSI/N dI/dt = βSI/N - γI dR/dt = γI ``` ## Biological Interpretation ### Disease Progression 1. **Susceptible (S)**: Individuals who can catch the disease 2. **Infected (I)**: Individuals who have the disease and can transmit it 3. **Recovered (R)**: Individuals who have recovered and gained immunity, or have been removed (died) ### Vaccination Effects - **Direct protection**: Vaccinated individuals move from S to R without getting infected - **Indirect protection**: Fewer susceptible individuals reduce transmission rate - **Herd immunity**: When enough people are immune, the disease cannot spread effectively ## Output Visualizations 1. **Time Series**: Shows S, I, R populations over time 2. **Phase Plane (S vs I)**: Trajectory showing relationship between susceptible and infected 3. **Peak Infection Comparison**: Effects of different initial conditions 4. **Vaccination Impact**: Comparison of scenarios with and without vaccination 5. **Cumulative Cases**: Total number of people who got infected over time 6. **3D Phase Space**: Complete trajectory through S-I-R space ## Key Insights ### Epidemic Dynamics - **Exponential growth phase**: Early rapid increase in infections - **Peak infections**: Maximum number of infected individuals - **Decline phase**: Infections decrease as susceptible population depletes - **Final epidemic size**: Total fraction of population that gets infected ### Vaccination Benefits - **Reduced peak infections**: Lower maximum number of simultaneous infections - **Delayed peak**: More time to prepare healthcare response - **Reduced total cases**: Fewer people get infected overall - **Flattened curve**: Spreads infections over longer time period ## Educational Applications This model helps students understand: - Compartmental models in epidemiology - Disease transmission dynamics - Public health intervention strategies - System of coupled differential equations - Phase space analysis - Parameter sensitivity analysis - Mathematical modeling of real-world phenomena - The importance of vaccination in disease control ## Codes ``` MATLAB % SIR Epidemic Model Driver % S' = -βSI/N - νS (Susceptible) % I' = βSI/N - γI (Infected) % R' = γI + νS (Recovered/Removed) % where N = S + I + R (total population) % Parameters beta = 0.5; % transmission rate gamma = 0.1; % recovery rate nu = 0.01; % vaccination rate N = 1000; % total population % Time span range = [0, 200]; h = 0.1; % Initial conditions [S, I, R] iv1 = [990, 10, 0]; % Initial outbreak scenario iv2 = [990, 20, 0]; % Larger initial outbreak % Define the SIR ODE system (row vector output) sir_ode = @(t, y) [-beta*y(1)*y(2)/N - nu*y(1), ... beta*y(1)*y(2)/N - gamma*y(2), ... gamma*y(2) + nu*y(1)]; % Solve using vectorRK4 for both scenarios [t1, y1] = vectorRK4(sir_ode, range, iv1, h); [t2, y2] = vectorRK4(sir_ode, range, iv2, h); % Also solve without vaccination for comparison sir_no_vac = @(t, y) [-beta*y(1)*y(2)/N, ... beta*y(1)*y(2)/N - gamma*y(2), ... gamma*y(2)]; [t3, y3] = vectorRK4(sir_no_vac, range, iv1, h); % Create plots figure; % Time series for first scenario subplot(2, 3, 1); plot(t1, y1(:, 1), 'b-', 'LineWidth', 2, 'DisplayName', 'Susceptible'); hold on; plot(t1, y1(:, 2), 'r-', 'LineWidth', 2, 'DisplayName', 'Infected'); plot(t1, y1(:, 3), 'g-', 'LineWidth', 2, 'DisplayName', 'Recovered'); title('SIR Model with Vaccination'); xlabel('Time (days)'); ylabel('Population'); legend; grid on; hold off; % Phase plane: S vs I subplot(2, 3, 2); plot(y1(:, 1), y1(:, 2), 'b-', 'LineWidth', 2, 'DisplayName', 'With Vaccination'); hold on; plot(y3(:, 1), y3(:, 2), 'r--', 'LineWidth', 2, 'DisplayName', 'No Vaccination'); plot(iv1(1), iv1(2), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b'); title('Phase Plane: S vs I'); xlabel('Susceptible'); ylabel('Infected'); legend; grid on; hold off; % Peak infection comparison subplot(2, 3, 3); plot(t1, y1(:, 2), 'b-', 'LineWidth', 2, 'DisplayName', 'I₀ = 10'); hold on; plot(t2, y2(:, 2), 'r--', 'LineWidth', 2, 'DisplayName', 'I₀ = 20'); title('Infected Population vs Initial Conditions'); xlabel('Time (days)'); ylabel('Infected Population'); legend; grid on; hold off; % Vaccination vs No vaccination subplot(2, 3, 4); plot(t1, y1(:, 2), 'b-', 'LineWidth', 2, 'DisplayName', 'With Vaccination'); hold on; plot(t3, y3(:, 2), 'r--', 'LineWidth', 2, 'DisplayName', 'No Vaccination'); title('Effect of Vaccination on Infections'); xlabel('Time (days)'); ylabel('Infected Population'); legend; grid on; hold off; % Cumulative cases (total who got infected) subplot(2, 3, 5); cumulative1 = N - y1(:, 1); % Total who left susceptible pool cumulative3 = N - y3(:, 1); plot(t1, cumulative1, 'b-', 'LineWidth', 2, 'DisplayName', 'With Vaccination'); hold on; plot(t3, cumulative3, 'r--', 'LineWidth', 2, 'DisplayName', 'No Vaccination'); title('Cumulative Cases'); xlabel('Time (days)'); ylabel('Total Cases'); legend; grid on; hold off; % 3D phase space (S, I, R) subplot(2, 3, 6); plot3(y1(:, 1), y1(:, 2), y1(:, 3), 'b-', 'LineWidth', 2, 'DisplayName', 'With Vaccination'); hold on; plot3(y3(:, 1), y3(:, 2), y3(:, 3), 'r--', 'LineWidth', 2, 'DisplayName', 'No Vaccination'); plot3(iv1(1), iv1(2), iv1(3), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b'); title('3D Phase Space (S, I, R)'); xlabel('Susceptible'); ylabel('Infected'); zlabel('Recovered'); legend; grid on; hold off; sgtitle('SIR Epidemic Model Analysis'); % Calculate and display key metrics R0 = beta / gamma; % Basic reproduction number peak_infected_vac = max(y1(:, 2)); peak_infected_no_vac = max(y3(:, 2)); final_recovered_vac = y1(end, 3); final_recovered_no_vac = y3(end, 3); fprintf('Basic Reproduction Number R₀ = %.2f\n', R0); fprintf('Peak Infected (with vaccination): %.0f people\n', peak_infected_vac); fprintf('Peak Infected (no vaccination): %.0f people\n', peak_infected_no_vac); fprintf('Final Recovered (with vaccination): %.0f people\n', final_recovered_vac); fprintf('Final Recovered (no vaccination): %.0f people\n', final_recovered_no_vac); fprintf('Lives saved by vaccination: %.0f people\n', final_recovered_no_vac - final_recovered_vac); ``` ``` Python """ SIR Epidemic Model with Vaccination Python implementation using scipy.integrate.solve_ivp and custom RK4/Euler methods """ import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp from mpl_toolkits.mplot3d import Axes3D def euler_sys(f, t_span, y0, h): """ Euler's method for systems of ODEs Parameters: f: function that takes (t, y) and returns dy/dt t_span: [t_start, t_end] y0: initial conditions array h: step size Returns: t: time array y: solution array (each column is one variable) """ t_start, t_end = t_span t = np.arange(t_start, t_end + h, h) n_steps = len(t) n_vars = len(y0) y = np.zeros((n_steps, n_vars)) y[0] = y0 for i in range(n_steps - 1): y[i+1] = y[i] + h * f(t[i], y[i]) return t, y def vector_rk4(f, t_span, y0, h): """ 4th-order Runge-Kutta method for systems of ODEs Parameters: f: function that takes (t, y) and returns dy/dt t_span: [t_start, t_end] y0: initial conditions array h: step size Returns: t: time array y: solution array (each column is one variable) """ t_start, t_end = t_span t = np.arange(t_start, t_end + h, h) n_steps = len(t) n_vars = len(y0) y = np.zeros((n_steps, n_vars)) y[0] = y0 for i in range(n_steps - 1): k1 = h * f(t[i], y[i]) k2 = h * f(t[i] + h/2, y[i] + k1/2) k3 = h * f(t[i] + h/2, y[i] + k2/2) k4 = h * f(t[i] + h, y[i] + k3) y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6 return t, y def main(): # Parameters beta = 0.5 # transmission rate gamma = 0.1 # recovery rate nu = 0.01 # vaccination rate N = 1000 # total population # Time span t_span = [0, 200] h = 0.1 # Initial conditions [S, I, R] iv1 = np.array([990, 10, 0]) # Initial outbreak scenario iv2 = np.array([990, 20, 0]) # Larger initial outbreak # Define the SIR ODE system with vaccination def sir_ode(t, y): """ SIR model with vaccination y[0] = S (Susceptible) y[1] = I (Infected) y[2] = R (Recovered) """ S, I, R = y dSdt = -beta * S * I / N - nu * S dIdt = beta * S * I / N - gamma * I dRdt = gamma * I + nu * S return np.array([dSdt, dIdt, dRdt]) # Define SIR system without vaccination for comparison def sir_no_vac(t, y): """ SIR model without vaccination """ S, I, R = y dSdt = -beta * S * I / N dIdt = beta * S * I / N - gamma * I dRdt = gamma * I return np.array([dSdt, dIdt, dRdt]) # Solve using custom RK4 for both scenarios t1, y1 = vector_rk4(sir_ode, t_span, iv1, h) t2, y2 = vector_rk4(sir_ode, t_span, iv2, h) # Solve without vaccination for comparison t3, y3 = vector_rk4(sir_no_vac, t_span, iv1, h) # Create plots fig = plt.figure(figsize=(15, 10)) # Time series for first scenario ax1 = plt.subplot(2, 3, 1) plt.plot(t1, y1[:, 0], 'b-', linewidth=2, label='Susceptible') plt.plot(t1, y1[:, 1], 'r-', linewidth=2, label='Infected') plt.plot(t1, y1[:, 2], 'g-', linewidth=2, label='Recovered') plt.title('SIR Model with Vaccination') plt.xlabel('Time (days)') plt.ylabel('Population') plt.legend() plt.grid(True) # Phase plane: S vs I ax2 = plt.subplot(2, 3, 2) plt.plot(y1[:, 0], y1[:, 1], 'b-', linewidth=2, label='With Vaccination') plt.plot(y3[:, 0], y3[:, 1], 'r--', linewidth=2, label='No Vaccination') plt.plot(iv1[0], iv1[1], 'bo', markersize=8, markerfacecolor='b') plt.title('Phase Plane: S vs I') plt.xlabel('Susceptible') plt.ylabel('Infected') plt.legend() plt.grid(True) # Peak infection comparison ax3 = plt.subplot(2, 3, 3) plt.plot(t1, y1[:, 1], 'b-', linewidth=2, label='I₀ = 10') plt.plot(t2, y2[:, 1], 'r--', linewidth=2, label='I₀ = 20') plt.title('Infected Population vs Initial Conditions') plt.xlabel('Time (days)') plt.ylabel('Infected Population') plt.legend() plt.grid(True) # Vaccination vs No vaccination ax4 = plt.subplot(2, 3, 4) plt.plot(t1, y1[:, 1], 'b-', linewidth=2, label='With Vaccination') plt.plot(t3, y3[:, 1], 'r--', linewidth=2, label='No Vaccination') plt.title('Effect of Vaccination on Infections') plt.xlabel('Time (days)') plt.ylabel('Infected Population') plt.legend() plt.grid(True) # Cumulative cases (total who got infected) ax5 = plt.subplot(2, 3, 5) cumulative1 = N - y1[:, 0] # Total who left susceptible pool cumulative3 = N - y3[:, 0] plt.plot(t1, cumulative1, 'b-', linewidth=2, label='With Vaccination') plt.plot(t3, cumulative3, 'r--', linewidth=2, label='No Vaccination') plt.title('Cumulative Cases') plt.xlabel('Time (days)') plt.ylabel('Total Cases') plt.legend() plt.grid(True) # 3D phase space (S, I, R) ax6 = plt.subplot(2, 3, 6, projection='3d') ax6.plot(y1[:, 0], y1[:, 1], y1[:, 2], 'b-', linewidth=2, label='With Vaccination') ax6.plot(y3[:, 0], y3[:, 1], y3[:, 2], 'r--', linewidth=2, label='No Vaccination') ax6.scatter([iv1[0]], [iv1[1]], [iv1[2]], color='b', s=50) ax6.set_title('3D Phase Space (S, I, R)') ax6.set_xlabel('Susceptible') ax6.set_ylabel('Infected') ax6.set_zlabel('Recovered') ax6.legend() plt.suptitle('SIR Epidemic Model Analysis') plt.tight_layout() plt.show() # Calculate and display key metrics R0 = beta / gamma # Basic reproduction number peak_infected_vac = np.max(y1[:, 1]) peak_infected_no_vac = np.max(y3[:, 1]) final_recovered_vac = y1[-1, 2] final_recovered_no_vac = y3[-1, 2] print(f'Basic Reproduction Number R₀ = {R0:.2f}') print(f'Peak Infected (with vaccination): {peak_infected_vac:.0f} people') print(f'Peak Infected (no vaccination): {peak_infected_no_vac:.0f} people') print(f'Final Recovered (with vaccination): {final_recovered_vac:.0f} people') print(f'Final Recovered (no vaccination): {final_recovered_no_vac:.0f} people') print(f'Lives saved by vaccination: {final_recovered_no_vac - final_recovered_vac:.0f} people') # Herd immunity threshold herd_immunity = 1 - 1/R0 print(f'Herd immunity threshold: {herd_immunity:.1%}') # Optional: Compare with scipy solve_ivp sol1 = solve_ivp(sir_ode, t_span, iv1, dense_output=True) t_scipy = np.linspace(t_span[0], t_span[1], 1000) y_scipy = sol1.sol(t_scipy).T plt.figure(figsize=(12, 8)) plt.subplot(2, 1, 1) plt.plot(t1, y1[:, 1], 'b-', linewidth=2, label='Custom RK4') plt.plot(t_scipy, y_scipy[:, 1], 'r--', alpha=0.7, label='SciPy solve_ivp') plt.title('Infected Population: Custom RK4 vs SciPy') plt.xlabel('Time (days)') plt.ylabel('Infected Population') plt.legend() plt.grid(True) plt.subplot(2, 1, 2) error = np.interp(t_scipy, t1, y1[:, 1]) - y_scipy[:, 1] plt.plot(t_scipy, error, 'g-', linewidth=1) plt.title('Difference between Custom RK4 and SciPy') plt.xlabel('Time (days)') plt.ylabel('Error in Infected Population') plt.grid(True) plt.tight_layout() plt.show() if __name__ == "__main__": main() ```