# Computational Models for Complex Systems Project

## Implementation of a SIRVSD model with different age goups

### The model
$$
\begin{equation}
    \begin{cases}
    \frac{\textit{dS}_j}{\textit{dt}} = \phi \textit{R}_j - \eta_j (t) \textit{S}_j + \rho \textit{V}_j - \textit{S}_j \sum_{k=1}^M \beta_{j,k} \textit{I}_k \\[4pt]
    \frac{\textit{dI}_j}{\textit{dt}} = - \gamma \textit{I}_j - \mu_j \textit{I}_j + \textit{S}_j \sum_{k=1}^M \beta_{j,k} \textit{I}_k \\[4pt]
    \frac{\textit{dR}_j}{\textit{dt}} = \gamma \textit{I}_j - \phi \textit{R}_j & \hspace{10pt} \forall j \in \left(0, M\right] \\[4pt]
    \frac{\textit{dV}_j}{\textit{dt}} = \eta_j (t) \textit{S}_j - \rho \textit{V}_j \\[4pt]
    \frac{\textit{dD}_j}{\textit{dt}} = \mu_j \textit{I}_j \\[4pt]
    \end{cases}
\end{equation}
$$

with $M$ = 4 age groups:
- Children (0-9 years)
- Teenagers (10-19 years)
- Adults (20-69 years)
- Senior Citizens (70+ years)

### Transition Schema


### Model assumptions
- Constant population and each compartment is normalized (S+I+R+V+D = 1).
- We don't consider vital dynamics (reproduction, "natural" death, migration), only the deaths caused by disease.
- Some coefficients are the same for all age group (only transmission, vaccination and mortality coefficients are different).
- We consider symmetric matrix for the transmission coefficients (heterogeneous contacts).
- Both recovery and vaccination immunity are not forever, but it ensures 100\% protection from infection.
- Just one injection of the vaccine is considered.
- Disease is transmitted by personal contacts between individuals of I and S classes (horizontal transmission).
- Contacts between individuals are random, the number of infections is proportional to both I and S.

### Definition of parameters
In our system, we have a lot of different coefficients:
- $\phi$ is the $\textit{transfer coefficient}$ for loss of immunity from Recovered.
- $\rho$ is the $\textit{transfer coefficient}$ for loss of immunity from Vaccinated.
- $\beta_{j, k}$ is the $\textit{infection coefficient}$, computed as the average number of contacts per person per time, multiplied by the probability of disease transmission in a contact between a susceptible subject of group $j$ and an infected subject of group $k$. <br>
We define also the entire $\textit{infection coefficient matrix}$ $\beta$ (or $\textit{contact matrix}$): 
$$\beta = \begin{bmatrix}
        \beta_{1, 1} & \cdots & \beta_{1, M}\\
        \vdots & \ddots & \vdots \\
        \beta_{M, 1} & \cdots & \beta_{M, M}
\end{bmatrix}
$$
- $\gamma$ is the $\textit{recovery coefficient}$ of each infected subject.
- $\mu_j$ is the $\textit{mortality coefficient}$, different for each age group.

The last coefficient is $\eta_j(t)$, a time-dependent $\textit{vaccination coefficient}$, defines as follows:
$$
    \begin{equation}
        \eta_j (t) = 
        \begin{cases}
            0 &\quad\text{if } t < t_{\textit{vacc}_j}\\
            \eta_j &\quad\text{otherwise}
        \end{cases}
    \end{equation}
$$

where $t_{\textit{vacc}_j}$ defines the starting day of the vaccination period.


```python
import numpy as np
from scipy.integrate import odeint, solve_ivp
import matplotlib.pyplot as plt
```


```python
def sirvd_solver(t,beta_matrix,gamma,mu_group,phi, rho,eta_group,x0,start_vaccination):
    """Wrapper function to compute ODEs using different APIs (currently only with methods of SciPy)

    Args:
        t (np.ndarray): simulation time
        beta_matrix (np.ndarray): infection coefficient for each group
        gamma (float): recovery coefficient (same for all group)
        mu_group (list): mortality coefficient for each group (case fatality rate ISS report January 2021)
        phi (float): transfer coefficient for loss of immunity from recovered (six months of immunity and same for all group)
        rho (float): transfer coefficient for loss of immunity from vaccinated (nine months of immunity and same for all group)
        eta_group (list): vaccination coefficient for each group
        x0 (list): initial conditions
        start_vaccination (list): day of start of the vaccination period for each group
    """
    def assign_vaccination_coefficient(t, eta, start_vaccination):
        """Auxiliary function to assign time-dependent vaccination coefficient eta

        Args:
            t (float): scalar representing the current timestamp
            eta (float): vaccination coefficient
            start_vaccination (int): starting day of vaccination

        Returns:
            float: eta for a specific timestamp (0 or eta)
        """
        if t < start_vaccination or start_vaccination == -1: # -1 means no vaccination for a specific age group
            return 0
        else:
            return eta
    
    def sirvd(t,x,beta_matrix,gamma,mu_group,phi,rho,eta_group,start_vaccination,assign_vaccination_coefficient):
        """
        Function called by solve_ivp (or odeint) to compute the derivative of x at t.
        """
        n_groups = len(start_vaccination) # or any other "_group" parameter
        derivatives_matrix = np.zeros((5,4)) # save all derivatives in a 2-D array
        n_infectious = [x[j+n_groups] for j in range(0,n_groups)] # list of the number of infectious measured for each group
        for j in range(0,n_groups):
            s = x[j] # Susceptible
            # i = n_infectious[j]
            i = x[j+n_groups] # Infectious
            r = x[j+n_groups*2] # Recovered
            v = x[j+n_groups*3] # Vaccinated
            d = x[j+n_groups*4] # Deceased
            eta = assign_vaccination_coefficient(t,eta_group[j],start_vaccination[j]) # time-dependent parameter
            derivatives_matrix[0][j] = phi*r - eta*s + rho*v - s*np.dot(beta_matrix[j],n_infectious) # dsdt
            derivatives_matrix[1][j] = s*np.dot(beta_matrix[j],n_infectious) - gamma*i - mu_group[j]*i # didt
            derivatives_matrix[2][j] = gamma*i - phi*r # drdt
            derivatives_matrix[3][j] = eta*s - rho*v # dvdt
            derivatives_matrix[4][j] = mu_group[j]*i # dddt
        return derivatives_matrix.reshape(-1) # return all measurements with a 1-D array

    # odeint solve a system of ordinary differential equations using lsoda from the FORTRAN library odepack.
    # y = odeint(sirvd,x0,t,tfirst=True,args=(beta_matrix,gamma,mu_group,phi,rho,eta_group,start_vaccination,assign_vaccination_coefficient,))
    # return y
    # for new code, use scipy.integrate.solve_ivp to solve a differential equation (SciPy documentation).
    sol = solve_ivp(sirvd,[t[0],t[-1]],x0,t_eval=t,args=(beta_matrix,gamma,mu_group,phi,rho,eta_group,start_vaccination,assign_vaccination_coefficient,))
    return sol.y.T
```

### Case Study: COVID-19 pandemic
### Definition of constants
We considered these groups:
- Children (0-9 years)
- Teenagers (10-19 years)
- Adults (20-69 years)
- Senior Citizens (70+ years)

In the following initial conditions and parameters, we define an indexing like this:

    S_0_GROUP = [x, y, w, y]
        - x is the intitial number of susceptible in "Children" group
        - y is the intitial number of susceptible in "Teenagers" group
        - w is the intitial number of susceptible in "Adults" group
        - z is the intitial number of susceptible in "Senior Citizens" group

We can change the value of the costants, according to the experiments we want to run and the plots we want to analyze.


```python
group_dict = {
    "children": 0,
    "teenagers": 1,
    "adults": 2,
    "senior": 3
}
START = 0 # observation starting day
END = 365 # observation ending day
START_VACCINATION_GROUP = [-1, -1, -1, -1] # day of start of the vaccination period (-1 means no vaccination)

# initial conditions
S_0_GROUP = [0.99, 0.99, 0.99, 0.99] # Susceptible
I_0_GROUP = [0.01, 0.01, 0.01, 0.01] # Infectious
R_0_GROUP = [0, 0, 0, 0] # Recovered
V_0_GROUP = [0, 0, 0, 0] # Vaccinated
D_0_GROUP = [0, 0, 0, 0] # Deceased

# model parameters
beta_matrix = np.array([[0.05,0.003,0.04,0.005],[0.003,0.09,0.07,0.007],[0.04,0.07,0.09,0.02],[0.005,0.007,0.02,0.03]]) # infection coefficient for each group
gamma = 1/15 # recovery coefficient (same for all group)
mu_group = [0.00009, 0.00005, 0.00688, 0.15987] # mortality coefficient for each group (case fatality rate ISS report January 2021)
phi = 1/180 # transfer coefficient for loss of immunity from recovered (six months of immunity and same for all group)
rho = 1/270 # transfer coefficient for loss of immunity from vaccinated (nine months of immunity and same for all group)
eta_group = [0.003, 0.003, 0.003, 0.01] # vaccination coefficient for each group

t = np.linspace(START,END,END-START+1) # setting the simulation time and the number of points
```

### Function call with same vaccination strategy


```python
results_dict = {}
x_0 = [*S_0_GROUP, *I_0_GROUP, *R_0_GROUP, *V_0_GROUP, *D_0_GROUP] # unpacking list operator
y = sirvd_solver(t, beta_matrix, gamma, mu_group, phi, rho, eta_group, x_0, START_VACCINATION_GROUP)
_, n_total_column = y.shape
n_groups = len(group_dict) # number of age groups
n_compartments = int(n_total_column/n_groups) # number of compartments of the model
for group_name, group_id in group_dict.items():
    # select the right columns (the compartments) for each age group 
    results_dict[group_name] = y[:,[group_id+n_groups*j for j in range(0,n_compartments)]]
```

### Function call with combination of different vaccination strategy


```python
"""vaccination_dict = {
        "no_vaccination": 0,
        "vaccination_strategy_ascending_order": 1,
        "vaccination_strategy_descending_order": 2,
        "vaccination_strategy_same_time": 3
    }
    results_dict = {}
    x_0 = [*S_0_GROUP, *I_0_GROUP, *R_0_GROUP, *V_0_GROUP, *D_0_GROUP] # unpacking list operator
    for vacc_name, vacc_id in vaccination_dict.items():
        results_dict[vacc_name] = {}
        if(vacc_id == 1): # ascending order
            START_VACCINATION_GROUP = [0, 30, 60, 90]
        elif(vacc_id == 2): # descending order
            START_VACCINATION_GROUP = [90, 60, 30, 0]
        elif(vacc_id == 3): # same time
            START_VACCINATION_GROUP = [0, 0, 0, 0]
            eta_group = [0.0025, 0.0025, 0.0025, 0.0025]
        y = sirvd_solver(t, beta_matrix, gamma, mu_group, phi, rho, eta_group, x_0, START_VACCINATION_GROUP)
        _, n_total_column = y.shape
        n_groups = len(group_dict) # number of age groups
        n_compartments = int(n_total_column/n_groups) # number of compartments of the model
        for group_name, group_id in group_dict.items():
            # select the right columns (the compartments) for each age group 
            results_dict[vacc_name][group_name] = y[:,[group_id+n_groups*j for j in range(0,n_compartments)]]"""
```

### References
- [Compartmental models in epidemiology — Wikipedia, the free encyclopedia](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology)
- [Modeling Infectious Diseases in Humans and Animals](https://homepages.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter3/Program_3.3/index.html)
- [Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models](https://link.springer.com/content/pdf/10.1007/BF00276956.pdf)
- [Analysis of COVID-19 Data with PRISM: Parameter Estimation and SIR Modelling](https://www.springerprofessional.de/en/analysis-of-covid-19-data-with-prism-parameter-estimation-and-si/18929902)
- [Use of a Modified SIRD Model to Analyze COVID-19 Data](https://pubs.acs.org/doi/pdf/10.1021/acs.iecr.0c04754)
- [Global  results  for  an  SIRS  model  with  vaccination and isolation](https://www.sciencedirect.com/science/article/abs/pii/S1468121810000763)
- [Mathematical  models  of  contact patterns between age groups for predicting the spread of infectious diseases](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4002176/pdf/nihms-570528.pdf)
- [A statistical methodology for data-driven partitioning of infectious disease incidence into age-groups](https://arxiv.org/pdf/1907.03441.pdf)
- [Lab24 - Coronavirus in Italia, i dati e la mappa](https://lab24.ilsole24ore.com/coronavirus/#)
- [SIR Modelling with data fitting in Python](https://github.com/Unipisa/SIR-covid)
- [Matplotlib Documentation](https://matplotlib.org/stable/contents.html)

495082f71cd76234f2a0118dd154738ac366dd65
ipynb
Jupyter Notebook
notebooks/Multi-age structured SIRVSD model.ipynb
Musca23/SIRVSD-model-with-age-groups
28a8587490ec99b0713b99d154ba4864a0517414
MIT
2021-12-22T16:47:13.000Z
2022-02-17T10:40:34.000Z
notebooks/Multi-age structured SIRVSD model.ipynb
Musca23/SIRVSD-model-with-age-groups
28a8587490ec99b0713b99d154ba4864a0517414
MIT
notebooks/Multi-age structured SIRVSD model.ipynb
Musca23/SIRVSD-model-with-age-groups
28a8587490ec99b0713b99d154ba4864a0517414
MIT
Qwen/Qwen-72B
1. YES
2. YES