Skip to content

Commit

Permalink
move to markdown
Browse files Browse the repository at this point in the history
  • Loading branch information
@anandijain @ChrisRackauckas
anandijain authored and ChrisRackauckas committed Jan 30, 2023
1 parent 8eb5555 commit c91c2c0
Show file tree
Hide file tree
Showing 2 changed files with 174 additions and 33 deletions.
18 changes: 7 additions & 11 deletions Scenario2/bin.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,16 +17,11 @@
# T, threatened (infected with life-threatening symptoms, detected);
# H, healed (recovered); E, extinct (dead).

cd(@__DIR__)
using OrdinaryDiffEq, ModelingToolkit, EasyModelAnalysis, SBML, SBMLToolkit, UnPack, Test
import Plots

@info "usings"

function sub_cont_ev(ev, rules)
ModelingToolkit.SymbolicContinuousCallback(substitute(ev.eqs, rules),
substitute(ev.affect, rules))
end
fn = "Scenario2/Giordano2020.xml"
fn = "Giordano2020.xml"

myread(fn) = readSBML(fn, doc -> begin
set_level_and_version(3, 2)(doc)
Expand Down Expand Up @@ -54,7 +49,6 @@ sys2 = ODESystem(eqs2, ModelingToolkit.get_iv(sys), states(sys), parameters(sys)
continuous_events = evs, defaults = defs, name = nameof(sys))
ssys = structural_simplify(sys2)
prob = ODEProblem(ssys, [], (0.0, 100.0))

sol = solve(prob, Tsit5())
@test sol.retcode == ReturnCode.Success
Plots.plot(sol)
Expand Down Expand Up @@ -98,6 +92,9 @@ xmax, xmaxval = get_max_t(prob_test1, sum(idart) * ITALY_POPULATION)
# Unit test #2:
# To reproduce Fig 2b and 2d, set Event_trigger_Fig3b = 0, Event_trigger_Fig3d = 0, Event_trigger_Fig4b = 0, Event_trigger_Fig4d = 0, epsilon_modifer = 1, alpha_modifer = 1 and run for t = 45 d
# -> These are already the settings in the Model


# V model
sysv = eval(quote
var"##iv#608" = (@variables(t))[1]
var"##sts#609" = (collect)(@variables(Infected(t), Healed(t), Extinct(t),
Expand Down Expand Up @@ -206,15 +203,14 @@ sysv = eval(quote
probv = ODEProblem(sysv, [], (0, 100))
solv = solve(probv, Tsit5())
plot(solv)
plot(solv, idxs=[og_states; Vaccinated])
plot(solv, idxs = [og_states; Vaccinated])
plot(solt1; idxs = sum(idart))

xmax, xmaxval = get_max_t(probv, sum(idart) * ITALY_POPULATION)
xmax, xmaxval = get_max_t(probv, sum(idart))

@test isapprox(xmax, 47; atol = 5)
@test isapprox(xmaxval, 0.6;atol=0.1)

@test isapprox(xmaxval, 0.6; atol = 0.1)

# @unpack Infected, Extinct = sys

Expand Down
189 changes: 167 additions & 22 deletions docs/src/Scenario2/Evaluation_Scenario_2.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,9 @@

```@example scenario2
cd(@__DIR__)
using OrdinaryDiffEq, ModelingToolkit, EasyModelAnalysis, SBML, SBMLToolkit, UnPack, Plots
using OrdinaryDiffEq, ModelingToolkit, EasyModelAnalysis, SBML, SBMLToolkit, UnPack, Test
import Plots
function sub_cont_ev(ev, rules)
ModelingToolkit.SymbolicContinuousCallback(substitute(ev.eqs, rules),
substitute(ev.affect, rules))
end
fn = "Giordano2020.xml"
myread(fn) = readSBML(fn, doc -> begin
Expand All @@ -27,25 +24,18 @@ defs_ = ModelingToolkit.defaults(sys)
defs = deepcopy(defs_)
evs = ModelingToolkit.get_continuous_events(sys)
# these are the constant=false params
params_to_sub = unique(ModelingToolkit.lhss(vcat(map(x -> x.affect, evs)...)))
@unpack alpha, epsilon, gamma, beta, delta, mu, nu, lambda, rho, kappa, xi, sigma, zeta, eta, theta = sys
@unpack Infected, Healed, Extinct, Diagnosed, Ailing, Recognized, Susceptible, Threatened = sys
@unpack alpha, epsilon, gamma, beta, delta, mu, nu, lambda, rho, kappa, xi, sigma, zeta, eta, theta, tau = sys
ps = [alpha, epsilon, gamma, beta, delta, mu, nu, lambda, rho, kappa, xi, sigma, zeta, eta]
@unpack Infected, Healed, Extinct, Diagnosed, Ailing, Recognized, Susceptible, Threatened = sys2
subd = Dict(params_to_sub .=> ps)
newevs = map(x -> sub_cont_ev(x, subd), evs)
@parameters t
D = Differential(t)
eqs2 = deepcopy(eqs)
append!(eqs2, D.(ps) .~ 0)
sys2 = ODESystem(eqs, ModelingToolkit.get_iv(sys), states(sys), parameters(sys);
continuous_events = newevs, defaults = defs, name = nameof(sys))
sys2 = structural_simplify(sys2)
```

```@example scenario2
prob = ODEProblem(sys2, [], (0.0, 1000.0))
sol = solve(prob, Tsit5())
plot(sol, idxs = Infected)
sys2 = ODESystem(eqs2, ModelingToolkit.get_iv(sys), states(sys), parameters(sys);
continuous_events = evs, defaults = defs, name = nameof(sys))
ssys = structural_simplify(sys2)
```

### Unit Tests
Expand All @@ -60,10 +50,50 @@ plot(sol, idxs = Infected)
> Simulate for 100 days, and determine the day and level of peak total infections (sum over all the infected states I, D, A, R, T). Expected output: The peak should occur around day 47, when ~60% of the population is infected.
```@example
```

#### Unit Test 2

> Now update the parameters to reflect various interventions that Italy implemented during the first wave, as described in detail on pg. 9. Simulate for 100 days, reproduce the trajectories in Fig. 2B, and determine the day and level of peak total infections (sum over all the infected states I, D, A, R, T). Expected output: Trajectories in Fig. 2B, peak occurs around day 50, with ~0.2% of the total population infected.
```@example
ITALY_POPULATION = 60e6
u0s = [
Infected => 200 / ITALY_POPULATION,
Diagnosed => 20 / ITALY_POPULATION,
Ailing => 1 / ITALY_POPULATION,
Recognized => 2 / ITALY_POPULATION,
Threatened => 0,
Healed => 0,
Extinct => 0,
]
push!(u0s, Susceptible => 1 - sum(last.(u0s)))
# The resulting basic reproduction number is R0 = 2.38.
pars = [alpha => 0.570, beta => 0.011, delta => 0.011, gamma => 0.456, epsilon => 0.171,
theta => 0.371,
zeta => 0.125, eta => 0.125, mu => 0.017, nu => 0.027, tau => 0.01,
lambda => 0.034, rho => 0.034, kappa => 0.017, xi => 0.017, sigma => 0.017]
prob = ODEProblem(ssys, [], (0, 100))
prob_test1 = remake(prob, u0 = u0s, p = pars)
solt1 = solve(prob_test1, Tsit5(); saveat = 0:100)
og_states = states(sys)[1:8]
idart = [Infected, Diagnosed, Ailing, Recognized, Threatened]
plot(solt1; idxs = Infected)
plot(solt1; idxs = Diagnosed)
plot(solt1; idxs = idart)
@test solt1[Infected + Healed] == solt1[Infected] + solt1[Healed]
plot(solt1.t, solt1[sum(idart)] * ITALY_POPULATION; label = "IDART absolute")
plot(solt1.t, solt1[sum(idart)]; label = "IDART percent")
xmax, xmaxval = get_max_t(prob_test1, sum(idart))
@test isapprox(xmax, 47; atol = 4)
@test isapprox(xmaxval, 0.002;atol=0.01)
```
### Sensitivity Analysis

> The difference between 1.b.i and 1.b.ii are changes in some parameter values over time. Describe the difference in outcomes between b.i and b.ii. Perform a sensitivity analysis to understand the sensitivity of the model to parameter variations and determine which parameter(s) were most responsible for the change in outcomes.
Expand Down Expand Up @@ -95,7 +125,122 @@ EasyModelAnalysis.optimal_parameter_intervention_for_threshold(prob, threshold_o
### Ingest SIDARTHE-V

```@example scenario2
sysv = eval(quote
var"##iv#608" = (@variables(t))[1]
var"##sts#609" = (collect)(@variables(Infected(t), Healed(t), Extinct(t),
Diagnosed(t), Ailing(t),
Recognized(t), Susceptible(t),
Threatened(t), Vaccinated(t),
alpha(t), epsilon(t), gamma(t),
beta(t), delta(t), mu(t), nu(t),
lambda(t), rho(t), kappa(t), xi(t),
sigma(t), zeta(t), eta(t)))
var"##ps#610" = (collect)(@parameters(ModelValue_21, epsilon_modifier, tau,
theta, ModelValue_19, ModelValue_20,
Event_trigger_Fig4b,
Event_trigger_Fig3d, ModelValue_16,
Event_trigger_Fig3b, alpha_modifier,
Event_trigger_Fig4d, ModelValue_17,
ModelValue_18, Italy, tau1, phi))
var"##eqs#611" = [(~)((Differential(t))(Infected),
(+)((*)((+)((/)((*)(Infected, alpha), Italy),
(/)((*)(Diagnosed, beta), Italy),
(/)((*)(Ailing, gamma), Italy),
(/)((*)(Recognized, delta), Italy)),
Susceptible), (*)(-1 // 1, Infected, epsilon),
(*)(-1 // 1, Infected, lambda),
(*)(-1 // 1, Infected, zeta)))
(~)((Differential(t))(Healed),
(+)((*)(Ailing, kappa), (*)(Diagnosed, rho),
(*)(Infected, lambda), (*)(Recognized, xi),
(*)(Threatened, sigma)))
(~)((Differential(t))(Extinct),
(*)(tau, Threatened) + tau1 * Recognized)
(~)((Differential(t))(Diagnosed),
(+)((*)(Infected, epsilon),
(*)(-1 // 1, Diagnosed, eta),
(*)(-1 // 1, Diagnosed, rho)))
(~)((Differential(t))(Ailing),
(+)((*)(Infected, zeta), (*)(-1 // 1, theta, Ailing),
(*)(-1 // 1, Ailing, kappa),
(*)(-1 // 1, Ailing, mu)))
(~)((Differential(t))(Recognized),
(+)((*)(theta, Ailing), (*)(Diagnosed, eta),
(*)(-1 // 1, Recognized, nu),
(*)(-1 // 1, Recognized, xi)))
(~)((Differential(t))(Susceptible),
(*)((+)((/)((*)(-1, Infected, alpha), Italy),
(/)((*)(-1, Diagnosed, beta), Italy),
(/)((*)(-1, Ailing, gamma), Italy),
(/)((*)(-1, Recognized, delta), Italy)) -
phi * Susceptible,
Susceptible))
(~)((Differential(t))(Threatened),
(+)((*)(Ailing, mu), (*)(Recognized, nu),
(*)(-1 // 1, tau, Threatened),
(*)(-1 // 1, Threatened, sigma)))
Differential(t)(Vaccinated) ~ phi * Susceptible
(~)((Differential(t))(alpha), -0.0);
(~)((Differential(t))(epsilon), -0.0);
(~)((Differential(t))(gamma), -0.0);
(~)((Differential(t))(beta), -0.0);
(~)((Differential(t))(delta), -0.0);
(~)((Differential(t))(mu), -0.0);
(~)((Differential(t))(nu), -0.0);
(~)((Differential(t))(lambda), -0.0);
(~)((Differential(t))(rho), -0.0);
(~)((Differential(t))(kappa), -0.0);
(~)((Differential(t))(xi), -0.0);
(~)((Differential(t))(sigma), -0.0);
(~)((Differential(t))(zeta), -0.0);
(~)((Differential(t))(eta), -0.0)]
var"##defs#612" = (Dict)((Pair)(delta, 0.011), (Pair)(xi, 0.017),
(Pair)(Diagnosed, 3.33333333e-7),
(Pair)(Event_trigger_Fig3b, 0.0),
(Pair)(Extinct, 0.0), (Pair)(kappa, 0.017),
(Pair)(zeta, 0.125), (Pair)(eta, 0.125),
(Pair)(nu, 0.027), (Pair)(Healed, 0.0),
(Pair)(Infected, 3.33333333e-6),
(Pair)(ModelValue_16, 0.0),
(Pair)(alpha_modifier, 1.0), (Pair)(Italy, 1.0),
(Pair)(Event_trigger_Fig3d, 0.0),
(Pair)(ModelValue_20, 1.0), (Pair)(sigma, 0.017),
(Pair)(Threatened, 0.0), (Pair)(lambda, 0.034),
(Pair)(alpha, 0.57),
(Pair)(Event_trigger_Fig4b, 0.0),
(Pair)(ModelValue_17, 0.0),
(Pair)(Event_trigger_Fig4d, 0.0),
(Pair)(Susceptible, 0.9999963),
(Pair)(beta, 0.011),
(Pair)(Recognized, 3.33333333e-8),
(Pair)(rho, 0.034), (Pair)(mu, 0.017),
(Pair)(epsilon, 0.171),
(Pair)(Ailing, 1.66666666e-8),
(Pair)(gamma, 0.456), (Pair)(ModelValue_19, 0.0),
(Pair)(ModelValue_21, 1.0), (Pair)(theta, 0.371),
(Pair)(epsilon_modifier, 1.0), (Pair)(tau, 0.01),
(Pair)(ModelValue_18, 0.0),
Vaccinated => 0,
tau1 => 0.0200,
phi => 0.0025)
var"##iv#613" = (@variables(t))[1]
(ODESystem)(var"##eqs#611", var"##iv#613", var"##sts#609", var"##ps#610";
defaults = var"##defs#612", name = Symbol("##SBML#530"),
checks = false)
end)
# todo set the event flags
# todo validate the new params
probv = ODEProblem(sysv, [], (0, 100))
solv = solve(probv, Tsit5())
plot(solv)
plot(solv, idxs = [og_states; Vaccinated])
plot(solt1; idxs = sum(idart))
xmax, xmaxval = get_max_t(probv, sum(idart) * ITALY_POPULATION)
xmax, xmaxval = get_max_t(probv, sum(idart))
@test isapprox(xmax, 47; atol = 5)
@test isapprox(xmaxval, 0.6; atol = 0.1)
```

### Setup the Parameters
Expand Down

0 comments on commit c91c2c0

Please sign in to comment.