Skip to content

Commit 6640e3b

Browse files
matthiaskoenigmatthiaskoenig
committed
work on sensitivity analysis
1 parent 4297ac7 commit 6640e3b

File tree

2 files changed

+187
-160
lines changed

2 files changed

+187
-160
lines changed

‎src/sbmlsim/sensitivity/global_sensitivity.py‎

Lines changed: 92 additions & 160 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,11 @@
11
"""Global sensitivity analysis.
22
33
TODO:
4-
- [ ] multi-output model
5-
- [ ] calculation of pharmacokinetic parameters (improved PK calculation)
6-
- [ ] get all parameter ids from model
4+
75
- [ ] get defined parameter bounds from model (annotate information);
8-
- [ ] storage of results
9-
- [ ] visualization of results
6+
- [ ] storage of results simulation
7+
- [ ] storage of results sensitivity analysis
8+
- [ ] visualization of results (heatmap)
109
- [ ] alternative methods:
1110
- [ ] local sensitivity analysis
1211
- [ ] SOBOL indices Sobol Sensitivity Analysis (Sobol 2001, Saltelli 2002, Saltelli et al. 2010)
@@ -31,17 +30,25 @@
3130
from pathlib import Path
3231
from dataclasses import dataclass
3332

34-
from pkdb_analysis.pk.pharmacokinetics import TimecoursePK
33+
3534
from sbmlutils.console import console
3635
from roadrunner._roadrunner import NamedArray
37-
from pint import UnitRegistry
36+
37+
3838

3939
@dataclass
40-
class SBMLSensitivityAnalysis:
41-
"""Parent class for sensitivity analysis."""
40+
class SensitivitySimulation:
41+
"""Base class for sensitivity calculation.
42+
43+
The sensitivity simulation runs a model simulation under a given set of
44+
model changes and returns a dictionary of scalar outputs.
45+
This function is called repeatedly during the sensitivity calculation.
46+
"""
47+
4248
model_path: Path
4349
selections: list[str]
4450
rr: roadrunner.RoadRunner = None
51+
outputs: list[str] = None
4552

4653
def __init__(self, model_path: Path, selections: list[str]):
4754
self.model_path = model_path
@@ -52,168 +59,93 @@ def __init__(self, model_path: Path, selections: list[str]):
5259
integrator.setSetting("variable_step_size", True)
5360
# state = rr.saveStateS()
5461

62+
# get the outputs from the simulation
63+
y = self.simulate(changes={})
64+
self.outputs = list(y.keys())
5565

56-
def losartan_simulation(self, changes: dict[str, float]) -> dict[str, float]:
57-
"""Run a given model simulation and create scalar readouts.
5866

59-
Applies the changes and gets the readouts. This is the general structure of the abstract
60-
run simulation. This can afterwards be parallelized.
67+
def simulate(self, changes: dict[str, float]) -> dict[str, float]:
68+
"""Runs a model simulation and returns the scalar results dictionary.
6169
62-
This must be a small subset of data
70+
This must be implemented by the subclass to work.
6371
"""
72+
raise NotImplemented
6473

65-
# rr = roadrunner.RoadRunner()
66-
# rr.loadStateS(state)
67-
self.rr.resetAll()
68-
all_changes = {
69-
"PODOSE_los": 10.0, # [mg]
70-
**changes
71-
}
72-
for key, value in all_changes.items():
74+
def plot(self) -> None:
75+
"""Plots the model simulation for debugging."""
76+
raise NotImplemented
77+
78+
def apply_changes(self, changes: dict[str, float], reset_all: bool=True) -> None:
79+
"""Apply changes after possible reset of the model."""
80+
if reset_all:
81+
self.rr.resetAll()
82+
for key, value in changes.items():
7383
# print(f"{key=} {value=}")
7484
self.rr.setValue(key, value)
7585

76-
s: NamedArray = self.rr.simulate(start=0, end=5 * 24 * 60)
77-
console.print(type(s))
78-
# self.plot_simulation(s)
79-
80-
y: dict[str, float] = self.losartan_pkpd_parameters(s)
81-
return y
82-
8386

8487

85-
def losartan_pkpd_parameters(self, s: NamedArray) -> dict[str, float]:
86-
"""This calculation is highly model dependent.
87-
88-
This requires units.
89-
"""
88+
@dataclass
89+
class SBMLSensitivityAnalysis:
90+
"""Parent class for sensitivity analysis."""
9091

91-
y: dict[str, float] = {}
92-
93-
# pharmacokinetics
94-
ureg = UnitRegistry()
95-
Q_ = ureg.Quantity
96-
97-
# losartan
98-
time = Q_(s["time"], "min")
99-
tcpk = TimecoursePK(
100-
time=time,
101-
concentration=Q_(s["[Cve_los]"], "mM"),
102-
substance="losartan",
103-
ureg=ureg,
104-
dose=Q_(10, "mg"),
105-
)
106-
pk_dict = tcpk.pk.to_dict()
107-
# console.print(pk_dict)
108-
for pk_key in [
109-
"aucinf",
110-
"cmax",
111-
"thalf",
112-
"kel",
113-
]:
114-
y[f"[Cve_los]_{pk_key}"] = pk_dict[pk_key]
115-
116-
# E3174, L158
117-
for sid in ["[Cve_e3174]", "[Cve_l158]"]:
118-
tcpk = TimecoursePK(
119-
time=time,
120-
concentration=Q_(s[sid], "mM"),
121-
substance="losartan",
122-
ureg=ureg,
123-
dose=None,
124-
)
125-
pk_dict = tcpk.pk.to_dict()
126-
# console.print(pk_dict)
127-
for pk_key in [
128-
"aucinf",
129-
"cmax",
130-
"thalf",
131-
"kel",
132-
]:
133-
y[f"{sid}_{pk_key}"] = pk_dict[pk_key]
134-
135-
# pharmacodynamics
136-
for sid, f in [
137-
("[ang1]", "max"),
138-
("[ang2]", "max"),
139-
("[ren]", "max"),
140-
("[ald]", "min"),
141-
("SBP", "min"),
142-
("DBP", "min"),
143-
("MAP", "min"),
144-
]:
145-
# minimal and maximal value of readout
146-
if f == "max":
147-
y[f"{sid}_max"] = np.min(s[sid])
148-
elif f == "min":
149-
y[f"{sid}_min"] = np.max(s[sid])
150-
151-
return y
152-
153-
def _plot_simulation(self, s: NamedArray) -> None:
154-
155-
# plotting
156-
from matplotlib import pyplot as plt
157-
plt.plot(
158-
s["time"], s["[Cve_los]"],
159-
# df["time"], df["[Cve_los]"],
160-
marker="o",
161-
markeredgecolor="black",
162-
color="tab:blue",
163-
)
164-
plt.show()
165-
166-
167-
def wrapped_run_simulation(self, X, func=losartan_simulation):
168-
# We transpose to obtain each column (the model factors) as separate variables
169-
changes: dict[str, float] = {}
170-
for k, key in enumerate(self.names):
171-
changes[key] = X[k]
172-
173-
# Then call the original model
174-
return list(func(self, changes).values())
175-
176-
177-
def calculate_sensitivity(self):
178-
179-
y = self.losartan_simulation(changes={})
180-
self.outputs = list(y.keys())
181-
self.names = ['BW']
182-
183-
# Defining the model inputs
184-
sp = ProblemSpec({
185-
'num_vars': len(self.names),
186-
'names': self.names,
187-
'bounds': [
188-
[50, 150],
189-
# [0.003, 0.005]
190-
],
191-
"outputs": self.outputs,
192-
})
193-
194-
# Generate samples
195-
samples = saltelli.sample(sp, 1024)
196-
sp.set_samples(samples)
197-
198-
199-
# Evaluate model
200-
# sp.evaluate(wrapped_run_simulation)
201-
202-
Y = np.zeros((samples.shape[0], len(self.outputs)))
203-
for k, X in enumerate(samples):
204-
print(k)
205-
Y[k, :] = self.wrapped_run_simulation(X)
206-
sp.set_results(Y)
207-
208-
209-
# Perform Analysis
210-
Si = sp.analyze(SALib.analyze.sobol)
211-
print(Si['S1'])
212-
print(Si['ST'])
213-
total_Si, first_Si, second_Si = Si.to_df()
214-
Si.plot()
215-
from matplotlib import pyplot as plt
216-
plt.show()
92+
sensitivity_simulation: SensitivitySimulation
93+
94+
def __init__(self, sensitivity_simulation: SensitivitySimulation):
95+
# assign simulation
96+
self.sensitivity_simulation = sensitivity_simulation
97+
98+
99+
# def wrapped_run_simulation(self, X, func=losartan_simulation):
100+
# # We transpose to obtain each column (the model factors) as separate variables
101+
# changes: dict[str, float] = {}
102+
# for k, key in enumerate(self.names):
103+
# changes[key] = X[k]
104+
#
105+
# # Then call the original model
106+
# return list(func(self, changes).values())
107+
108+
109+
# def calculate_sensitivity(self):
110+
#
111+
# y = self.losartan_simulation(changes={})
112+
# self.outputs = list(y.keys())
113+
# self.names = ['BW']
114+
#
115+
# # Defining the model inputs
116+
# sp = ProblemSpec({
117+
# 'num_vars': len(self.names),
118+
# 'names': self.names,
119+
# 'bounds': [
120+
# [50, 150],
121+
# # [0.003, 0.005]
122+
# ],
123+
# "outputs": self.outputs,
124+
# })
125+
#
126+
# # Generate samples
127+
# samples = saltelli.sample(sp, 1024)
128+
# sp.set_samples(samples)
129+
#
130+
#
131+
# # Evaluate model
132+
# # sp.evaluate(wrapped_run_simulation)
133+
#
134+
# Y = np.zeros((samples.shape[0], len(self.outputs)))
135+
# for k, X in enumerate(samples):
136+
# print(k)
137+
# Y[k, :] = self.wrapped_run_simulation(X)
138+
# sp.set_results(Y)
139+
#
140+
#
141+
# # Perform Analysis
142+
# Si = sp.analyze(SALib.analyze.sobol)
143+
# print(Si['S1'])
144+
# print(Si['ST'])
145+
# total_Si, first_Si, second_Si = Si.to_df()
146+
# Si.plot()
147+
# from matplotlib import pyplot as plt
148+
# plt.show()
217149

218150

219151

‎src/sbmlsim/sensitivity/parameters.py‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
"""Tools and helpers to handle parameters for senitivity analysis."""
2+
from pathlib import Path
3+
from typing import Optional
4+
5+
import libsbml
6+
import numpy as np
7+
from sbmlutils.console import console
8+
9+
def bounds_for_parameters():
10+
"""Retrieve the parameter bounds from the model."""
11+
12+
13+
14+
def parameters_for_sensitivity_analysis(
15+
sbml_path: Path,
16+
exclude_ids: Optional[set[str]] = None,
17+
exclude_na: bool = True,
18+
) -> dict[str, str]:
19+
"""Retrieve parameters from model for the sensitivity analysis.
20+
21+
Constant parameters, constant compartments and constant species are returned.
22+
23+
:sbml_path: Path to the SBML file.
24+
:param exclude_ids: ids to exclude,
25+
:param exclude_na: whether to exclude NA values
26+
:return: dict[id, name]
27+
"""
28+
29+
30+
doc: libsbml.SBMLDocument = libsbml.readSBMLFromFile(str(sbml_path))
31+
sbml_model: libsbml.Model = doc.getModel()
32+
33+
id2name: dict[str, str] = {}
34+
excluded: set[str] = set()
35+
36+
# constant parameters
37+
p: libsbml.Parameter
38+
for p in sbml_model.getListOfParameters():
39+
sid = p.getId()
40+
if p.getConstant() is True:
41+
if exclude_na and np.isnan(p.getValue()):
42+
excluded.add(sid)
43+
continue
44+
id2name[sid] = p.getName() if p.isSetName() else sid
45+
46+
# constant compartments
47+
c: libsbml.Compartment
48+
for c in sbml_model.getListOfCompartments():
49+
sid = c.getId()
50+
if c.getConstant() is True:
51+
if exclude_na and np.isnan(c.getSize()):
52+
excluded.add(sid)
53+
continue
54+
id2name[sid] = c.getName() if c.isSetName() else sid
55+
56+
# constant species or boundaryCondition == True
57+
s: libsbml.Species
58+
for s in sbml_model.getListOfSpecies():
59+
sid = s.getId()
60+
name = s.getName() if s.isSetName() else sid
61+
if exclude_na:
62+
if not s.isSetInitialAmount() and not s.isSetInitialConcentration():
63+
excluded.add(sid)
64+
continue
65+
if s.isSetInitialAmount() and np.isnan(s.getInitialAmount()):
66+
excluded.add(sid)
67+
continue
68+
if s.isSetInitialConcentration() and np.isnan(s.getInitialConcentration()):
69+
excluded.add(sid)
70+
continue
71+
72+
if s.getConstant() is True or s.getBoundaryCondition() is True:
73+
id2name[sid] = name
74+
75+
# remove excluded ids
76+
if exclude_ids:
77+
for sid in exclude_ids:
78+
if sid in id2name:
79+
excluded.add(sid)
80+
id2name.pop(sid)
81+
82+
console.print(f"Excluded parameters: {excluded}")
83+
84+
return id2name
85+
86+
if __name__ == "__main__":
87+
model_path = Path(__file__).parent / "models" / "losartan" / "losartan_body_flat.xml"
88+
id2name = parameters_for_sensitivity_analysis(
89+
sbml_path=model_path,
90+
exclude_ids = {
91+
"conversion_min_per_day", # constant conversion factor
92+
"Mr_los", # molecular weight
93+
}
94+
)
95+
console.print(id2name)

0 commit comments

Comments
 (0)