Skip to content

Commit

Permalink
Add initial implementation of scenario 2
Browse files Browse the repository at this point in the history
  • Loading branch information
@bgyori
bgyori committed Jul 7, 2023
1 parent e24e626 commit e1822b8
Show file tree
Hide file tree
Showing 2 changed files with 1,679 additions and 0 deletions.
286 changes: 286 additions & 0 deletions notebooks/hackathon_2023.07/scenario2.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,286 @@
{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"id": "7798c424",
"metadata": {},
"source": [
"# Hackathon Scenario 2: Multiple Vaccines and Variants\n",
"Models from early in the Covid-19 pandemic were relatively simple, as we had no vaccines available, only NPIs to mitigate the spread of the disease, and only one strain (the original wild type) to worry about. By mid-2021, the situation was much more complex, as we had access to multiple vaccines, each with multiple doses, and the emergence of multiple variants. You are interested in modeling this situation with SV2(AIR)3 (https://www.nature.com/articles/s41598- 022-06159-x, https://doi.org/10.1038/s41598-022-06159-x), one of the first models to represent multiple Covid-19 variants, and the competition between them."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "d274fc77",
"metadata": {},
"source": [
"1. Extract the simplest version of the model with 1 vaccine and 1 variant (shown in Fig. 1A)\n",
"\n",
" a. (TA1 Model Extraction Workflow) Extract the model from Matlab code\n",
"‘svair_simple.m’ and ‘get_beta.m’ and with the publication equations and Fig. 1A (with the exception of one small difference between code and figure). Note that in this code, the units of all state variables are in terms of fraction of the total population.\n",
" \n",
" b. (TA3 Simulation Workflow) Simulate the first 14 months of the pandemic beginning on January 1st, 2020. Consider only the Wild Type variant, and no vaccination during this period. Use parameter values from the supplementary material, for the wild type. Use the population of Ontario in 2019: 14.57 million people, and let the natural death rate 𝜇 = 2.05𝑒-5. For initial conditions, let 𝐼(0) = 1e-6, 𝑆(0) = 1 − 𝑠𝑢𝑚(𝐼(0) + 𝐴(0)), and all other values be 0. When do I(t) and A(t) peak, and what is the value of these variables at their peaks? How do the I(t) and A(t) profiles compare with Fig. 3d,f (which considered 3 variants and 2 types of vaccines)? Given that this question assumes only the presence of the Wild Type variant, do the results seem reasonable?\n",
" \n",
" c. (TA2 Domain Knowledge Grounding; ASKEM Workbench Only) Demonstrate that the domain knowledge groundings accurately reflect new concepts with the extensions such as “variants”."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "0991209a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'\\nS = y(ind); ind=ind+1; % original susceptible\\nSVR = y(ind); ind=ind+1; % lost immunity after vaccination or recovery\\nV1 = y(ind:ind+n_vax-1); ind=ind+n_vax; % one-dose vaccination\\nV2 = y(ind:ind+n_vax-1); ind=ind+n_vax; % fully vaccinated\\nI = y(ind:ind+n_var-1); ind=ind+n_var; % infected\\nIV = y(ind:ind+n_var-1); ind=ind+n_var; % infected even with vaccination\\nIR = y(ind:ind+n_var-1); ind=ind+n_var; % infected again after recovery from a different variant\\nA = y(ind:ind+n_var-1); ind=ind+n_var; % asymptomatic infections\\nAR = y(ind:ind+n_var-1); ind=ind+n_var; % asymptomatic infections after recovery from a different variant\\nR = y(ind:ind+n_var-1); ind=ind+n_var; % recovered\\nR2 = y(ind); ind=ind+1; % recovered after getting both variants\\n\\n\\n'"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import sympy\n",
"from mira.metamodel import *\n",
"from mira.modeling import Model\n",
"from mira.modeling.askenet.petrinet import AskeNetPetriNetModel\n",
"\n",
"\"\"\"\n",
"S = y(ind); ind=ind+1; % original susceptible\n",
"SVR = y(ind); ind=ind+1; % lost immunity after vaccination or recovery\n",
"V1 = y(ind:ind+n_vax-1); ind=ind+n_vax; % one-dose vaccination\n",
"V2 = y(ind:ind+n_vax-1); ind=ind+n_vax; % fully vaccinated\n",
"I = y(ind:ind+n_var-1); ind=ind+n_var; % infected\n",
"IV = y(ind:ind+n_var-1); ind=ind+n_var; % infected even with vaccination\n",
"IR = y(ind:ind+n_var-1); ind=ind+n_var; % infected again after recovery from a different variant\n",
"A = y(ind:ind+n_var-1); ind=ind+n_var; % asymptomatic infections\n",
"AR = y(ind:ind+n_var-1); ind=ind+n_var; % asymptomatic infections after recovery from a different variant\n",
"R = y(ind:ind+n_var-1); ind=ind+n_var; % recovered\n",
"R2 = y(ind); ind=ind+1; % recovered after getting both variants\n",
"\n",
"\n",
"\"\"\""
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "716217a7",
"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",
"populations = ['S', 'SVR', 'V1', 'V2', 'I', 'IV', 'IR',\n",
" 'A', 'AR', 'R', 'R2']\n",
"\n",
"c = {pop: Concept(name=pop, units=dimensionless_units)\n",
" for pop in populations}\n",
"\n",
"parameters = [\n",
" Parameter(name='beta', value=3.3e-9),\n",
" Parameter(name='beta_v1', value=0.2*3.3e-9),\n",
" Parameter(name='beta_v2', value=0.05*3.3e-9),\n",
" Parameter(name='beta_R', value=0.05*3.3e-9),\n",
" Parameter(name='gamma', value=1/28),\n",
" Parameter(name='mu', value=298.7),\n",
" Parameter(name='mu_I', value=0.001),\n",
" Parameter(name='mu_IV', value=0.15*0.001),\n",
" Parameter(name='nu_R', value=1/(4*365)),\n",
" Parameter(name='nu_v1', value=1/365),\n",
" Parameter(name='nu_v2', value=1/(4*365)),\n",
" Parameter(name='ai', value=0.5), # paper uses sigma instead of ai\n",
" Parameter(name='ai_V', value=0.85),\n",
" Parameter(name='ai_R', value=0.85),\n",
" Parameter(name='ai_beta_ratio', value=3.0)\n",
"]\n",
"\n",
"initials = {\n",
" 'S': Initial(concept=c['S'], value=1-1e-6),\n",
" 'SVR': Initial(concept=c['SVR'], value=0),\n",
" 'V1': Initial(concept=c['V1'], value=0),\n",
" 'V2': Initial(concept=c['V2'], value=0),\n",
" 'I': Initial(concept=c['I'], value=1e-6),\n",
" 'IV': Initial(concept=c['IV'], value=0),\n",
" 'IR': Initial(concept=c['IR'], value=0),\n",
" 'A': Initial(concept=c['A'], value=0),\n",
" 'AR': Initial(concept=c['AR'], value=0),\n",
" 'R': Initial(concept=c['R'], value=0),\n",
" 'R2': Initial(concept=c['R2'], value=0),\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "f91e5e06",
"metadata": {},
"outputs": [],
"source": [
"S, SVR, V1, V2, I, IV, IR, A, AR, R, R2 = sympy.symbols('S SVR V1 V2 I IV IR A AR R R2')"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "45045c21",
"metadata": {},
"outputs": [],
"source": [
"beta, beta_v1, beta_v2, beta_R, gamma, mu, mu_I, mu_IV, nu_R, nu_v1, nu_v2, ai, ai_V, ai_R, ai_beta_ratio = \\\n",
" sympy.symbols('beta beta_v1 beta_v2 beta_R gamma mu mu_I mu_IV nu_R nu_v1 nu_v2 ai ai_V ai_R ai_beta_ratio')"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "563bef26",
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
"dSdt = \n",
"- sum(beta.*S.*(I_total)) \n",
"- sum(vr1.*S) \n",
"+ mu*(1-S);\n",
"\n",
"dSVRdt = \n",
"+ nu_v1.*sum(V1) \n",
"+ nu_v2.*sum(V2) \n",
"+ sum(nu_R.*R) \n",
"+ nu_R.*R2 \n",
"- sum(beta.*SVR.*(I_total)) \n",
"- mu*SVR;\n",
"\n",
"dV1dt = \n",
"+ vr1.*(S+sum(A)) \n",
"- vr2.*V1 \n",
"- nu_v1.*V1 \n",
"- sum(beta_v1.*(V1*(I_total)'),2) \n",
"- mu*V1;\n",
"\n",
"dV2dt = + vr2.*V1 - nu_v2.*V2 - sum(beta_v2.*(V2.*(I_total)'),2) - mu*V2;\n",
"dIdt = (1-ai).*(+ beta.*S.*(I_total) + beta.*SVR.*(I_total)) - gamma.*I - mu_I.*I;\n",
"dIVdt = (1-ai_V).*(+ sum(beta_v1.*(V1*(I_total)'))' + sum(beta_v2.*(V2.*(I_total)'))') - gamma.*IV - mu_IV.*IV;\n",
"dIRdt = (1-ai_R).*(+ beta_R.*Rv.*(I_total)) - gamma.*IR - mu*IR;\n",
"dAdt = ai.*(+ beta.*S.*(I_total) + beta.*SVR.*(I_total)) + ai_V.*(sum(beta_v1.*(V1*(I_total)'))' + sum(beta_v2.*(V2.*(I_total)'))') - sum(vr1).*A - gamma.*A - mu.*A;\n",
"dARdt = ai_R.*(+ beta_R.*Rv.*(I_total)) - gamma.*AR - mu*AR;\n",
"dRdt = + gamma.*(I_total) - nu_R.*R - beta_R.*Rv.*(I_total) - mu*R;\n",
"dR2dt = + sum(gamma.*(IR+AR)) - nu_R.*R2 - mu*R2;\n",
"\n",
"\"\"\"\n",
"\n",
"all_infectious = c['I'], c['IV'], c['IR'], c['A'], c['AR']\n",
"infection_rate_term = I + IV + IR + ai_beta_ratio * (A + AR)\n",
"\n",
"templates = [\n",
" # Infection asym\n",
" GroupedControlledConversion(subject=c['S'], outcome=c['A'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta*S*infection_rate_term),\n",
" GroupedControlledConversion(subject=c['V1'], outcome=c['A'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta_v1*V1*infection_rate_term),\n",
" GroupedControlledConversion(subject=c['V2'], outcome=c['A'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta_v2*V2*infection_rate_term),\n",
"\n",
" # Infection sym\n",
" GroupedControlledConversion(subject=c['S'], outcome=c['I'],\n",
" controllers=all_infectious,\n",
" rate_law=(1-ai)*beta*S*infection_rate_term),\n",
" GroupedControlledConversion(subject=c['V1'], outcome=c['IV'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta_v1*V1*infection_rate_term),\n",
" GroupedControlledConversion(subject=c['V2'], outcome=c['IV'],\n",
" controllers=all_infectious,\n",
" rate_law=ai*beta_v2*V2*infection_rate_term),\n",
" # Vaccination\n",
" # Note: structurally we want to have these in the model but we make the rates 0\n",
" NaturalConversion(subject=c['S'], outcome=c['V1'], rate_law=0*S),\n",
" NaturalConversion(subject=c['A'], outcome=c['V1'], rate_law=0*A),\n",
" NaturalConversion(subject=c['V1'], outcome=c['V2'], rate_law=0*V1),\n",
" # Recovery\n",
" NaturalConversion(subject=c['I'], outcome=c['R'], rate_law=gamma*I),\n",
" NaturalConversion(subject=c['IV'], outcome=c['R'], rate_law=gamma*IV),\n",
" NaturalConversion(subject=c['IR'], outcome=c['R'], rate_law=gamma*IR),\n",
" NaturalConversion(subject=c['A'], outcome=c['R'], rate_law=gamma*A),\n",
" NaturalConversion(subject=c['AR'], outcome=c['R'], rate_law=gamma*AR),\n",
"\n",
" # Reinfection asym\n",
" GroupedControlledConversion(subject=c['R'], outcome=c['AR'],\n",
" controllers=all_infectious,\n",
" rate_law=ai_R*beta_R*R*infection_rate_term),\n",
" # Reinfection sym\n",
" GroupedControlledConversion(subject=c['R'], outcome=c['IR'],\n",
" controllers=all_infectious,\n",
" rate_law=(1-ai_R)*beta_R*R*infection_rate_term),\n",
"\n",
" # Second recovery\n",
" NaturalConversion(subject=c['IR'], outcome=c['R2'], rate_law=gamma*IR),\n",
" NaturalConversion(subject=c['AR'], outcome=c['R2'], rate_law=gamma*AR),\n",
"\n",
" # Death\n",
" NaturalDegradation(subject=c['S'], rate_law=mu*S),\n",
" NaturalDegradation(subject=c['SVR'], rate_law=mu*SVR),\n",
" NaturalDegradation(subject=c['V1'], rate_law=mu*V1),\n",
" NaturalDegradation(subject=c['V2'], rate_law=mu*V2),\n",
" NaturalDegradation(subject=c['I'], rate_law=mu_I*I),\n",
" NaturalDegradation(subject=c['IV'], rate_law=mu_IV*IV),\n",
" NaturalDegradation(subject=c['IR'], rate_law=mu*IR),\n",
" NaturalDegradation(subject=c['A'], rate_law=mu*A),\n",
" NaturalDegradation(subject=c['AR'], rate_law=mu*AR),\n",
" NaturalDegradation(subject=c['R'], rate_law=mu*R),\n",
" NaturalDegradation(subject=c['R2'], rate_law=mu*R2),\n",
"\n",
" # Loss of immunity\n",
" NaturalConversion(subject=c['V1'], outcome=c['SVR'], rate_law=nu_v1*V1),\n",
" NaturalConversion(subject=c['V2'], outcome=c['SVR'], rate_law=nu_v2*V2),\n",
" NaturalConversion(subject=c['R'], outcome=c['SVR'], rate_law=nu_R*R),\n",
" NaturalConversion(subject=c['R2'], outcome=c['SVR'], rate_law=nu_R*R2),\n",
"]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "25dfad3d",
"metadata": {},
"outputs": [],
"source": [
"tm = TemplateModel(templates=templates\n",
" parameters={p.name: p for p in parameters},\n",
" initial_conditions=initials,\n",
"AskeNetPetriNetModel(Model(tm)).to_json_file('scenario2_a.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 e1822b8

Please sign in to comment.