SAMENVATTING ModSim |Practicum

Modelleren en Simuleren Biosystemen|PRACTICA SAMENVATTING
Gegeven door Prof. Michiel Stock


------------ Ordinary Differential Equations(ODEs) --------------
!!in een wiskundig model gebruik je ~ om een relatie te beschrijven ipv =
▪ Variabelen definiëren @species
▪ Constanten definiëren @parameters → vast
▪ Parameters definiëren @parameters → vast of laten veranderen (leeg)

Netwerk → Problem → Solve → Plot


Modelopbouw: Catalyst
using ModelingToolKit: t_nounits as t, D_nounits as D
Model_name = @reaction_network begin
@species S(t) I(t) thr = waarde |niet perse invullen bhl bij threshold!!
@parameters 𝛼 𝛽 |niet perse invullen, constanten wel

#de reactie: Snelheid (boven pijl), [substraat]→ [product]
𝛼 ∗ 𝛽, S+I --> 2I
r*m, I --> D
end


 @species
o De t geeft weer van wat de waarde afhankelijk is (tijd hier)
 Kan de ODE’s tonen via convert(ODESystem, Model_name)
 Kan ook <--> gebruiken bv (k1, k1), 2NO <--> N2O2 |notatie voor beide kanten
 𝛼 ∗ 𝛽 komt boven de pijl en duid op de snelheid (senlehidsparameter k)



Van netwerk naar simulatie
Een netwerk is slechts een beschrijving, nog omzetten naar een ODEProblem

Stap 1: ingrediënten klaarzetten

u0 = [:S => waarde, :I => waarde, …] of [y=>1, D(y) =>0] #startwaarden
tspan = (0.0, 90.0) #over welke tijd
params = [:𝛼 => waarde, :𝛽 => waarde, :r => r]


 In het initiële model kennen we bij @params dus geen waarden toe, dat kunnen we
later doen
 Notatie :X wil zeggen dat Julia de ‘X’ als symbool moet interpreteren
 r willen we laten variëren dus zullen we een slider toevoegen waarover we r
kunnen laten variëren via de macro @bind


@bind r Slider(0.1:0.1:1, default=0.1, show_value=true)

 Variëert van 0.1 → 1 in stappen van 0.1 met startwaarde 0.1 en toont getal

,Stap 2: probleem oplossen = simuleren

oprob = ODEProblem(model_naam, u0, tspan, params) |probleem aanmaken wiskundig
osol = solve(oprob, Tsit5(), saveat = 0.5) |probleem oplossen

 Tsit5() is het algoritme dat de berekening uitvoert voor een ODE
 osol bevat nu alle data van de simulatie
 saveat = 0.5 → zegt op welke tijdstippen (0.5s) je de resultaten (alle
variabelen) wil opslaan in het osol object dus bij tspan = (0,10) krijg je 21
datapunten → het zal dus de grafiek bepalen


Condities en Events
Soms veranderd het systeem (parameters tijdens de simulatie (bv een lockdown op dag
14) → om dit toe te voegen gebruik je discrete events

▪ Tijds-event (discrete_events!)
o conditie1 =[14.0] => [𝛽 ~ 𝛽/2] #vn dag 14 halveert parameter 𝛽
▪ Treshold (continuous_events!)
o conditie2 =[S ~ thr] => [S ~ S – 500, thr ~ 1e9] #zdr S thr bereikt S-500
o Wordt in @species gedefinieerd in het model
o In een treshold conditie zorgen dat je de treshold waarde opnieuw aanpast
zodat je dit event niet meer zal meemaken (eenmalig)

!!Belangrijk: toevoegen event aan systeem
Als je een event toevoegt, moet je het model afmaken via complete() voordat je het
kunt oplossen

@named model_c = ReactionSystem(equations(model_naam), discrete_events=condition1)

o equations(model_naam) pakt reactieregels (S + I -> 2I) uit eerder gemaakt model
o discrete_events=condition1 plakt de extra regels aan het model

Model_complete = complete(model_c) #geüpdate model met conditie erin verwerkt

 Macro @named om het systeem een identiteit te geven
 Hierna opnieuw o_prob_c = ODEProblem(Model_complete, u0, tspan, params) …
 Bij events moet je deepcopy toevoegen om oorspronkelijk systeem te behouden


Handige hulp
@unpack S, I = model_naam

 Kan in het vervolg gewoon S en I gebruiken ipv telkens model_name.S
zelfde voor parameters en de thr waarde!
 Indien niet moest je in paragraaf hierboven model_name.S gebruiken en
model_name.𝛽 en model_name.thr

deepcopy(o_prob)

 Systeem kan soms aanpassen als je een simulatie runt, als je events/condities
toevoegt, via dit bouw je een geheugen in zodat je een exacte kopie hebt van
het originele model voor een volgende simulatie

oprob_c = ODEproblem(Model_complete, u0, tspan, params)
osol_c = solve(deepcopy(oprob_c), Tsit5(), saveat=0.5); #zodat het niet opslaat

, Steady state
Steady state = toestand waarin systeem niet meer veranderd (aanvoer = afvoer)
Kijken dus niet meer naar een tspan maar naar het punt waar alle afgeleiden = 0

 SteadyStateProblem
 u_guess

We gebruiken de eindwaarden van de simulatie om in u_guess in te vullen want grote
kans dat die al dichter bij het evenwichtpunt zullen liggen (osol[:S][end] bv)

u_guess = [:S => osol[:S][end], :X => osol[:X][end]]
sol_eq = solve(SteadyStateProblem(model_naam, u_guess, params))



Plotten
Gewone plots

plot(osol) |tekent alle variabelen (S, I, …) over tijd
plot(osol, idxs =[:S, :I]) |tekent enkel variabele S en I met t op x-as
plot(osol, idxs =(:S, :I)) |S op x-as, I op y-as ronde haken! → relatie tss
osol[:S][end] |geeft de laatste waarde van variabele S op einde

 Argumenten bij plot()
o osol #resultaten simulatie
o idxs =[ ] of ( ) #welke variabelen gebruiken
o xlab = “naam”, ylab = “naam” #asnamen
o ylim =(0, 100) of xlim = (0, 200)
o lw = 0.5 |line width

Meerdere simulaties in één figuur → plot!()

 ! in julia betekend pas het bestaande aan

begin

plot(osol1)
plot!(osol2) |voegt data van osol2 toe aan de grafiek die er al staat

end

 Argumenten bij plot!()
o linestyle = :dash |maakt stippellijn van de lijn om te onderscheiden
o label = none |zorgt dat toegevoegde lijn niet in legende verschijnt
o color = :grey |lijnkleur

Histogram: frequenties (aantal voorkomen bepaalde waarde)

histogram(times, bins=range(0, 120, length=121), xlab="Age (years)",
ylab="Frequency")

 Argumenten bij histogram() (ook histogram!() mogelijk)
o times = inputdata
o bins =range(begin_x, eind_x, length = aantal bins)
o color = :blue
o alpha = 0.5 |maakt staven half transparant (2 over elkaar bv)
o label = “Groep A” |geeft een naam toe voor de legend

, Modelopbouw: MTK algemeen
Catalyst (in MTK) → @reaction_network → voor chemische reacties
ModelingToolKit → lost algemene natuurkundige en wiskundige vergelijkingen op
𝑑
▪ D = afgeleide naar de tijd (𝑑𝑡)
▪ D(y) = y’ (snelheid)
▪ D(D(y) = y” (versnelling)




Verschillen

 Ipv een pijl (--> of bv <-->) gebruik je een tilde ~
o eq = m*D(D(y)) ~ -m*g -k*(y-l) – nu * D(y) |merk: F = m*a = …
o eq = 𝑚 ∗ 𝑦" ~ − 𝑚𝑔 − 𝑘(𝑦 − 𝑙) − 𝜈 ∗ 𝑦′
 Gebruikt dus voor het opstellen van een model dat niet een “reactie” is
 @variables y(t) ipv @species
 We stapelen nu de onderdelen op elkaar tot een systeem
 Catalyst waar alles binnen de begin en end van het netwerk staat


using ModelingToolKit: t_nounits as t, D_nounits as D

1. #Definieer onafhankelijke variabele (tijd)
@variables t
D = Differential(t) #of in using D_nounits as D

2. #Definieer veranderlijke variabelen en de parameters
@variables y(t)
@parameters m=2.0 g=9.81 k l=1.0 nu=0.5 |hoeven niet perse allemaal ingevuld

3. #Definieer de vergelijking (de logica)
eq = m * D(D(y)) ~ -m*g - k*(y - l) - nu*D(y)

4. #Bouw het systeem
@mtkbuild model = ODESystem(eq, t)


---- rest zelfde als catalyst ----

u0 = [ y => 1.5, # Startpositie D(y) => 0.0 # Startsnelheid] #initiële condities
tspan = (0.0, 20.0)

5. #Probleem oplossen
oprob = ODEProblem(model, u0, tspan) # Parameters zitten al in de @parameters regel
#je kan wel nog ( …, [k => 1.2, nu => 0.2]) toevoegen (parameters overschrijven)
osol = solve(oprob, Tsit5(), saveat = 0.5)

6. #Visualisatie
plot(osol, idxs=[y])

, Examples Hoofdstuk 2
Example 5: afkoelen koffie + melk erbij (p36)
Afkoelen koffie in ruimte op kamertemperatuur waarbij na 2 seconden melk wordt
bijgegoten (event)

using ModelingToolKit; T_nounits as t, D_nounits as D

@variables V(t) T(t)
@parameters Cp rho Te Tm k q_in

volume_balance = D(V) ~ q_in
heat_balance = expand_deriviatives(D(V*T)) ~ q_in*Tm + k*V*(Te – T)

---------------------------------
Expand_deriviatives is al een advanced gebruik van MTK
Gebruikt om symbolisch de afgeleide van een product van twee variabelen uit te werken
Schrijft automatisch V*D(T) + T*D(V)
---------------------------------

eqs = [volume_balance, heat_balance] #geeft de twee vgl hierin

dV, dT = symbolic_lineair_solve(eqs, [D(V), D(T)]
---------------------------------
Door gebruik van expand_derivitiatives vouwt Julia de gecombineerde energiebalans
open, hierna is een symbolische solver: symbolic_lineair_solve instaat om de losse
afgeleides D(V) en D(T) te isoleren en het ODESystem correct op te bouwen
---------------------------------

#Event maken van melk gieten

melk_gieten = [
[2] => [q_in ~ 5e-2] #op tijd 2 krijgt q_in een waarde (start gieten)
[3] => [q_in ~ 0] #op tijd 3 krijgt q_in de waarde 0 (stop gieten)
]


@mtkbuild kof_model = ODEsystem([D(V)~dV, D(T)~dT], t, discrete_events = melk_gieten)

u0 = [V=>1.5e-1, T=>80],
tspan = (0.0, 8.0)
#parameters invullen
kof_prob = ODEProblem(kof_model, u0, tspan, [Te => 20, Tm => 5, k => 0.1, q_in => 0]
kof_sol = solve(kof_prob)

plot(kof_sol)



Example 13: Populatie van 2 species met competitie (logistische groei) (p53)
Species A met max groeisnel. ra en species B met rb en draagkracht K en mortaliteit m

two_species_model = @reaction_network begin
@parameters ra = 0.2, rb = 0.15, K = 100, m = 0.1
ra*(1-(A+B)/K), A --> 2A
rb*(1-(A+B)/K), B --> 2B
m, (A,B) --> 0
end

,Example 12: Repressie en gen expressie (Reverse Hill vergelijking) (p51)
Repressor-concentratie (R/L) stijgt waardoor de gen-expressie afneemt (v)
Wat dus wil zeggen dat de hoeveelheid mRNA afneemt (hier eigen expressie in vb)

repressor_model = @reaction_network begin
@species R(t) = 0 mRNA(t) = 0

hillr(R, vtranscr, Ki, 4), 0 --> mRNA #transcriptie reactie
vtransl, mRNA --> mRNA + R #translatie reactie
d, mRNA --> 0 #degradatie mRNA
r, R --> 0 #degradatie van R
end


Compartiment modellen → complex systeem naar onderling verbonden compartimenten

Example 14: Farmacokinetiek (geneesmiddelen in het lichaam)
Geneesmiddel in systeem lichaam → ingenomen in maag (G) + geabsorbeerd in bloed (B)
+ komt terecht in weefsel (T)

pharkin = @reaction_network begin
@species B(t)=0 T(t)=0 #Startconcentraties in bloed en weefsel
k1, G --> B
k2, B --> T
end


 Pijlen stellen geen chemische reactie voor maar fysiek transport naar elkaar
 k1 en k2, de snelheden, bepalen de uitwisseling (eerste orde processen)
 Startwaarde van G later in ODEProblem definiëren


Example 15: Celcompartimenten (Gen-expressie)
Systeem eukaryotische cel → compartimenten: celkern en cytosol
Transcriptie in cel van gen → translatie in cytoplasma met vorming dimeer-eiwitten
(regulator dat transcriptie kan blokken)
#Transcriptie en regulatie in nucleus
nuc = @network_component nuc begin
alpha, G --> G + M
(p/V,m), D + G <--> DG
end

#Translatie en dimerisatie in cytosol
cyto = @network_component cyto begin

beta, M --> M + P
(k/V,k), 2P <--> D
sigma, P --> 0
mu, M --> 0
end

#Difinieer transport tussen beide
cell_model = @network_component model begin
gamma, $(nuc.M) --> $(cyto.M)
delta, $(cyto.D) --> $(nuc.D)
end

@named cell_model = compose(cell_model, [nuc, cyto] #hangt de drie submodellen aan elkaar

,Example 16: Reactor met dode zone
Bij aannemen dat chemische reactor niet perfect gemengd is
→ niet voldaan aan homogeniteitsassumptie
Opvangen door op te splitsen in goed gemengde bulkzone (Am) en dode zone (Ad)
De dode zone vormt een fractie f van het totaal volume f = Ad/(Am + Ad)
Wel uitwisseling tussen beide → beschrijven

reactor2 = @reaction_network begin
@species Ad(t)=0 B(t)=0
@parameters V=100 k=0.3 r=0.05 q=1 c=.1

q*c, 0 --> Am # Instroom in Am → nulde-orde
q/(V*(1-f)), Am --> 0 # Uitstroom in Am, gecorrigeerd door volume

r, Am --> B # Reactie die in het vat optreed (enkel in Am)

(k/(V*(1-f)), k/(V*f)), Am <--> Ad # uitwisseling met de dode zone
end


 Reactiesnelheden zijn hier niet gewoon meer parameters maar volledige vgl


Example 17: SIR-model (epidemologie met gezondheidscomparimenten)

sir = @reaction_network begin
@species S(t)=50 I(t)=5 R(t)=0
beta, S + I --> 2I
gamma, I --> R
end



Example 18: Vlinder-model (ecologie met levencycluscomparimenten)
Een vlinder kent verschillende stadia: ei(E), rups(C), kokon(P) en vlinder
Met de vlinder die een nieuw e ikan leggen om de cyclus te herstarten
Maar elke levensfase heeft een kans op sterven
Overgang (s), sterven (m), vruchtbaarheid (f)

butterfly = @reaction_network begin
@species E(t)=10 C(t)=0 P(t)=0 B(t)=0
@parameters f=1.5 s=0.05 m=0.05

s, E --> C
s, C --> P
s, P --> B

m, (E, C, P, B) --> 0 #manier #reacties samen schrijven

f, B --> B + E
end

, Algemene opmerkingen

• Bij het kiezen van welke afhankelijke variabele je wil onderzoeken, kan je bv
in een chemische reactie voor concentratie kiezen (mol/L) of voor Volume
terwijl een concentratie afhangt van het volume dus moeilijk in een behoudswet
te steken is → je kan ervoor kiezen beide in @variables/@species te steken, het
model dat gegenereerd wordt zal wanneer mogelijk de overbodige variabelen
verwijderen
• Symbolics pakket
o Gebruikt voor symbolische wiskunde waarmee je afgeleiden kan berekenen




o expand_derivatives() voert de afgeleiden uit op de achtergrond
o eval(build_function(df_sym, x)) om in Juliafunctie om te zetten

• Automatic differentiation: soort computer-gegenereerde afgeleiden

Forward mode (ForwardDiff) Reverse mode (Zygote of Enzyme)

,------------ Stochastic Differential Equations(SDEs) ------------
→ er zit randomness in model
SDE = model voor continuevariabelen, waar een expliciete ruiscomponent is toegevoegd
De vergelijking bestaat nu uit een deterministische drift (ODE) + een willekeurige
diffusie (Wiener proces gebouwd op Brownse beweging)

 ODE Tsit5() |kijkt naar gladde vergelijking
 SDE SRIW1() |houdt rekening met extra term voor ruis (Brownse beweging)
EM()


Monte Carlo simulatie is de methode die je hierop toepast:

1. Je voert de simulatie met SRIW1 bijvoorbeeld 1000 keer uit.
2. Elke run is uniek door de willekeurige ruis
3. Je berekent het gemiddelde en de standaardafwijking van al die 1000 runs.

Dit proces noemen we in Julia een Ensemble Problem. Je gebruikt Monte Carlo om de
statistiek te begrijpen van de stochastische vergelijkingen die SRIW1 voor je oplost.



Modelbouw SDE (continue ruis)
Model_naam = @reaction_network begin
@parameters η=40 #ruis
@default_noise_scaling η #stelt algemene ruis in op parameter η
α * β, S + I --> 2I, [noise_scaling = 60.0] #voor deze specifieke reactie wil
je een Andere ruis → overschrijft η
r * m, I --> D #gebruikt default ruis
r * (1 - m), I --> R
end

 Aan elke reactive/vergelijking wordt nu ruis η toegevoegd, dit kan aangepast
worden via [noise_scaling = waarde] voor een specifieke reactive

 Als je [noise_scaling = 0.0] zet, maak je die reactive puur deterministisch
(zonder willekeurg)

 Hebben nu ons model opgebouwd met een stelsel SDEs



Netwerk → simuleren
sprob = SDEProblem(model_naam, u0, tspan, params) #ipv ODEProblem
ssol = solve(sprob, EM() of SRIW1(), dt=0.1)

 Hebben nu conceptueel network omgezet naar wiskundig probleem dat kan opgelost
worden

 Hebben nu de waarden van één run, maar elke run is uniek en willikeurig pad +
geeft een andere grafiek (oplossen via ensemble simulatie)

 Maar! Omdat ruis in een grafiek wild kan uitslaan en willekeurig is, kan een
waarde in Julia onder 0 dalen (biologisch niet mogelijk) of veel te hoog worden
Oplossen via DiscreteCallback

, end Callbacks bij instabiliteit

Definiëren wanneer ingegrepen moet worden
function condition(u, t, integrator)

true #dwingt Julia na elke tijdsstap te kijken of er iets mis is

end

 u stelt “huidige staat” voor (een array met de waarden voor S, I, D en R)
 t stelt “huidige tijd” voor in de simulatie
 integrator is een object dat alles bijhoudt
o integrator.u |integrator.t |integrator.p = huidige parameters)
 Werken bv niet met t == 5 want weten niet op welk moment het onder 0 zal duiken
of te hoog worden want random → daarom via true (na elke stap checken)


Dwingen van de integratie om de waarden tussen 0 – 10.000.000 te houden
Vertellen wat er moet veranderen als de conditie bereikt is

function affect!(integrator)
for i = 1:4 #want hier vier variabelen (S,I,D en R) in u
if integrator.u[i] > 10000000 #max
integrator.u[i] = 10000000
end
if integrator.u[i] < 0 #min
integrator.u[i] = 0
end
end
end

 De ! achter een functie in Julia geeft aan dat de invoer overschreven moet
worden, dus als over 1000000 wordt gegaan → waarde op 100000 stellen (en 0)


Samenvoegen van de twee functies tot één regel om aan de solver te voeren
cb = DiscreteCallback(condition, affect!, save_positions=(false,false))

 save_positions=(false,false) zegt dat de correctie stil op achtegrond
uitgevoerd moeten worden zonder extra datapunten op te slaan op de grafiek
om werkgeheugen niet vol te laten lopen



Run → ensemble simulatie (Monte Carlo)
Nu ensembleProblem maken (bv om 100x simulatie te doorlopen en variate te bekijken)
en oplossen met callback in verwerkt
esprob = EnsembleProblem(sprob)
essol = solve(esprob, EM(), dt=0.1, callback=cb, trajectories=100)

 Kan nog argument save_everystep = true aanvullen
o Dwingt pc om exact elke tijdsstap (0.1) data op te slaan (vloeiend)
o Want soms slaat solver enkel begin, eind en sommige data ertussen op
 trajectories vertelt om het pad 100x volledig onafh te simuleren