Skip to content

Commit

Permalink
Add scenario 1 implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
@bgyori
bgyori committed Jun 30, 2023
1 parent 79a11a9 commit cb223c7
Show file tree
Hide file tree
Showing 4 changed files with 1,002 additions and 0 deletions.
244 changes: 244 additions & 0 deletions notebooks/hackathon_2023.07/scenario1.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,244 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "7a646367",
"metadata": {},
"source": [
"### Step a-b: create the base model\n",
"Questions: You want to implement a masking intervention in a simple compartmental model and simulate epidemic trajectories under different compliance scenarios. You found an existing model that incorporates masking as a time-dependent modification to the β parameter (https://doi.org/10.3390/ijerph18179027 plus accompanying code), and you want to ensure that the model is working as expected by reproducing plots in the publication. \n",
"\n",
"Replicate an analysis from the paper.\n",
"- (TA1 Model Extraction Workflow, TA2 Model Representation) Extract the SEIRD model (equations 1-5, with time-varying β as defined in equations 8-9) and load into the workbench. In equation 8, let kappa = γR0, R0 = 5. In equation 9, let k = 5.\n",
"\n",
"- (TA3 Simulation Workflow, Unit Test): Replicate Figure 3 from the paper, which maps to the first scenario in the paper- implementing a masking intervention at several different timepoints (delays of 0 days, 50 days, 100 days, and control case, from the date of first infection) in the pandemic, with 100% compliance. Recreate the Fig. 3 curves (including peak infection times and levels), up to some reasonable margin of error. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "56172b31",
"metadata": {},
"outputs": [],
"source": [
"import sympy\n",
"from mira.metamodel import *\n",
"from mira.modeling import Model\n",
"from mira.modeling.askenet.petrinet import AskeNetPetriNetModel"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "f07383fb",
"metadata": {},
"outputs": [],
"source": [
"person_units = Unit(expression=sympy.Symbol('person'))\n",
"day_units = Unit(expression=sympy.Symbol('day'))\n",
"per_day_units = Unit(expression=1/sympy.Symbol('day'))\n",
"dimensionless_units = Unit(expression=sympy.Integer('1'))\n",
"\n",
"c = {\n",
" 'S': Concept(name='S', units=dimensionless_units),\n",
" 'E': Concept(name='E', units=dimensionless_units),\n",
" 'I': Concept(name='I', units=dimensionless_units),\n",
" 'R': Concept(name='R', units=dimensionless_units),\n",
" 'D': Concept(name='D', units=dimensionless_units)\n",
"}\n",
"\n",
"\n",
"parameters = {\n",
" 'gamma': Parameter(name='gamma', value=1/11, units=per_day_units),\n",
" 'delta': Parameter(name='delta', value=1/5, units=per_day_units),\n",
" 'alpha': Parameter(name='alpha', value=0.000064, units=dimensionless_units),\n",
" 'rho': Parameter(name='rho', value=1/9, units=per_day_units),\n",
" 'N': Parameter(name='N', value=5_600_000, units=person_units),\n",
" 'beta_s': Parameter(name='beta_s', value=1),\n",
" 'beta_c': Parameter(name='beta_c', value=0.4),\n",
" 't_0': Parameter(name='t_0', value=89),\n",
" # D=11, gamma = 1/D, R_0 = 5 and\n",
" # beta = R_0 * gamma * mask(t) so kappa = 5/11\n",
" 'kappa': Parameter(name='kappa', value=5/11),\n",
"}\n",
"\n",
"initials = {\n",
" 'S': Initial(concept=Concept(name='S'), value=5_600_000-1),\n",
" 'E': Initial(concept=Concept(name='E'), value=1),\n",
" 'I': Initial(concept=Concept(name='I'), value=0),\n",
" 'R': Initial(concept=Concept(name='R'), value=0),\n",
" 'D': Initial(concept=Concept(name='D'), value=0),\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "bdf70b6a",
"metadata": {},
"outputs": [],
"source": [
"S, E, I, R, D, N, kappa, beta_s, beta_c, k, t_0, t, alpha, delta, rho, gamma = \\\n",
" sympy.symbols('S E I R D N kappa beta_s beta_c k t_0 t alpha delta rho gamma')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "d1b732b3",
"metadata": {},
"outputs": [],
"source": [
"m_1 = (beta_s - beta_c) / (1 + sympy.exp(-k*(t_0-t))) + beta_c\n",
"beta = kappa*m_1\n",
"\n",
"t1 = ControlledConversion(subject=c['S'],\n",
" outcome=c['E'],\n",
" controller=c['E'],\n",
" rate_law=S*I*beta / N)\n",
"t2 = NaturalConversion(subject=c['E'],\n",
" outcome=c['I'],\n",
" rate_law=delta*E)\n",
"t3 = NaturalConversion(subject=c['I'],\n",
" outcome=c['R'],\n",
" rate_law=(1-alpha)*gamma*I)\n",
"t4 = NaturalConversion(subject=c['I'],\n",
" outcome=c['D'],\n",
" rate_law=alpha*rho*I)\n",
"templates = [t1, t2, t3, t4]\n",
"tm = TemplateModel(\n",
" templates=templates,\n",
" parameters=parameters,\n",
" initials=initials,\n",
" time=Time(name='t', units=day_units)\n",
")\n",
"AskeNetPetriNetModel(Model(tm)).to_json_file('scenario1_a.json')"
]
},
{
"cell_type": "markdown",
"id": "4c29d51a",
"metadata": {},
"source": [
"### Step c: update beta\n",
"\n",
"(TA2 Model Modification Workflow, TA3 Simulation Workflow): Update the β(t) function to be defined as equations 8 and 10, with k1 = 5, and k2 = 1. This reflects the paper’s second scenario, gradual noncompliance with the masking policy over time. Rerun the simulation with several different delays in enforcing a mask policy (ranging from 0 to 140 days), and replicate Fig. 5, up to some reasonable margin of error."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "3413db92",
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"k_1, k_2, t_1, beta_nc = sympy.symbols('k_1 k_2 t_1 beta_nc')\n",
"parameters['k_1'] = Parameter(name='k_1')\n",
"parameters['k_2'] = Parameter(name='k_2')\n",
"parameters['t_1'] = Parameter(name='t_1', value=154)\n",
"parameters['beta_nc'] = Parameter(name='beta_nc', value=0.5)\n",
"m_2 = (beta_s - beta_c) / (1 + sympy.exp(-k_1*(t_0-t))) + (beta_c - beta_nc) / (1 + sympy.exp(-k_2*(t_1-t))) + beta_nc"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "f487b47f",
"metadata": {},
"outputs": [],
"source": [
"t1.rate_law = SympyExprStr(kappa*m_2)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "30a47395",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\kappa \\left(\\beta_{nc} + \\frac{\\beta_{c} - \\beta_{nc}}{1 + e^{- k_{2} \\left(- t + t_{1}\\right)}} + \\frac{- \\beta_{c} + \\beta_{s}}{1 + e^{- k_{1} \\left(- t + t_{0}\\right)}}\\right)$"
],
"text/plain": [
"kappa*(beta_nc + (beta_c - beta_nc)/(1 + exp(-k_2*(-t + t_1))) + (-beta_c + beta_s)/(1 + exp(-k_1*(-t + t_0))))"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t1.rate_law.args[0]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "0791ae2d",
"metadata": {},
"outputs": [],
"source": [
"AskeNetPetriNetModel(Model(tm)).to_json_file('scenario1_c.json')"
]
},
{
"cell_type": "markdown",
"id": "5f67f2e9",
"metadata": {},
"source": [
"### Step d: update for reinfection\n",
"(TA2 Model Modification Workflow, TA3 Simulation Workflow): Update the system of equations to include equations 6 and 7. This adds the potential for reinfection. Compare with the outcomes from 1c. What impact does immunity loss and potential for reinfection have?"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "38b6002e",
"metadata": {},
"outputs": [],
"source": [
"epsilon = sympy.Symbol('epsilon')\n",
"parameters['epsilon'] = Parameter(name='epsilon', value=1/90)\n",
"t5 = NaturalConversion(subject=c['R'],\n",
" outcome=c['S'],\n",
" rate_law=epsilon*R)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "3aed70ad",
"metadata": {},
"outputs": [],
"source": [
"AskeNetPetriNetModel(Model(tm)).to_json_file('scenario1_d.json')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Loading

0 comments on commit cb223c7

Please sign in to comment.