Skip to content

Commit

Permalink
Refactor volumes and heart rate
Browse files Browse the repository at this point in the history
  • Loading branch information
@finsberg
finsberg committed Oct 12, 2024
1 parent b3acaf3 commit 3b46667
Show file tree
Hide file tree
Showing 5 changed files with 49 additions and 44 deletions.
30 changes: 19 additions & 11 deletions examples/regazzoni_bleeding.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
setup_logging()



# Run first Zenker to get the correct heart rate for normal conditions
zenker_normal = Zenker()
zenker_normal.solve(T=100.0, dt=1e-3, dt_eval=0.1)
Expand All @@ -23,7 +24,7 @@


# Now we will simulate a bleeding and compute a new heart rate
blood_loss_parameters = {"start_withdrawal": 1, "end_withdrawal": 2, "flow_withdrawal": -2000, "flow_infusion": 0}
blood_loss_parameters = {"start_withdrawal": 1, "end_withdrawal": 200, "flow_withdrawal": -100, "flow_infusion": 0}
zenker_bleed = Zenker(parameters=blood_loss_parameters)
zenker_bleed.solve(T=300.0, dt=1e-3, dt_eval=0.1, initial_state=zenker_normal.state)
HR_bleed = zenker_bleed.results["fHR"][-1]
Expand Down Expand Up @@ -57,7 +58,7 @@
regazzoni_bleed_parmeters["circulation"]["external"] = blood_loss_parameters

regazzoni_bleed = Regazzoni2020(parameters=regazzoni_bleed_parmeters)
regazzoni_bleed.solve(num_cycles=20, initial_state=regazzoni_normal.state, dt_eval=dt_eval)
regazzoni_bleed.solve(num_cycles=50, initial_state=regazzoni_normal.state, dt_eval=dt_eval)
regazzoni_bleed.print_info()
N_bleed = int(regazzoni_bleed.HR / dt_eval)

Expand All @@ -75,18 +76,25 @@


fig, ax = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(10, 5))
ax[0, 0].plot(regazzoni_normal.results["V_LV"], regazzoni_normal.results["p_LV"])
ax[0, 0].set_xlabel("V [mL]")
ax[0, 0].set_ylabel("p [mmHg]")
ax[0, 0].set_title("Normal")
ax[1, 0].plot(regazzoni_normal.results["V_LV"][-N_normal:], regazzoni_normal.results["p_LV"][-N_normal:])
ax[1, 0].set_xlabel("V [mL]")
ax[0, 1].plot(regazzoni_bleed.results["V_LV"], regazzoni_bleed.results["p_LV"])
ax[0, 1].set_xlabel("V [mL]")
ax[0, 1].set_ylabel("p [mmHg]")
ax[0, 0].plot(regazzoni_normal.results["V_LV"], regazzoni_normal.results["p_LV"], label="LV")
ax[0, 0].plot(regazzoni_normal.results["V_RV"], regazzoni_normal.results["p_RV"], label="RV")
ax[1, 0].plot(regazzoni_normal.results["V_LV"][-N_normal:], regazzoni_normal.results["p_LV"][-N_normal:], label="LV")
ax[1, 0].plot(regazzoni_normal.results["V_RV"][-N_normal:], regazzoni_normal.results["p_RV"][-N_normal:], label="RV")
ax[0, 1].set_title("Bleeding")
ax[1, 1].plot(regazzoni_bleed.results["V_LV"][-N_bleed:], regazzoni_bleed.results["p_LV"][-N_bleed:])
ax[0, 1].plot(regazzoni_bleed.results["V_LV"], regazzoni_bleed.results["p_LV"], label="LV")
ax[0, 1].plot(regazzoni_bleed.results["V_RV"], regazzoni_bleed.results["p_RV"], label="RV")
ax[1, 1].plot(regazzoni_bleed.results["V_LV"][-N_bleed:], regazzoni_bleed.results["p_LV"][-N_bleed:], label="LV")
ax[1, 1].plot(regazzoni_bleed.results["V_RV"][-N_bleed:], regazzoni_bleed.results["p_RV"][-N_bleed:], label="RV")

ax[0, 0].set_ylabel("p [mmHg]")
ax[1, 0].set_ylabel("p [mmHg]")
ax[1, 0].set_xlabel("V [mL]")
ax[1, 1].set_xlabel("V [mL]")

for axi in ax.flatten():
axi.grid()
axi.legend()
plt.show()


Expand Down
13 changes: 5 additions & 8 deletions src/circulation/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,7 @@ def time_varying_elastance(self, EA, EB, tC, TC, TR, **kwargs):
tC=tC,
TC=TC,
TR=TR,
HR=self.HR,
RR=1 / self.HR,
)

def flux_through_valve(self, p1: float, p2: float, R: Callable[[float, float], float]) -> float:
Expand Down Expand Up @@ -233,7 +233,7 @@ def solve(
logger.info("Running circulation model")
if T is None:
assert num_cycles is not None, "Please provide num_cycles or T"
T = self.HR * num_cycles
T = num_cycles / self.HR

initial_state = initial_state or dict()

Expand All @@ -257,24 +257,21 @@ def solve(
else:
self.state = remove_units(self.state)

self.store(t)

time_start = time.time()

i = 0
while t < T:
self.callback(self, t, i % output_every_n_steps == 0)
self.step(t, dt)
if i % output_every_n_steps == 0:
self.store(t)
self.callback(self, t, i % output_every_n_steps == 0)

self.step(t, dt)
if self._verbose:
self.print_info()

if i % checkoint_every_n_steps == 0:
self.save_state()
t += dt
i += 1

duration = time.time() - time_start

logger.info(f"Done running circulation model in {duration:.2f} s")
Expand Down
29 changes: 16 additions & 13 deletions src/circulation/regazzoni2020.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -296,27 +296,30 @@ def step(self, t, dt):

@property
def volumes(self):
C_VEN_SYS = self.parameters["circulation"]["SYS"]["C_VEN"]
C_AR_SYS = self.parameters["circulation"]["SYS"]["C_AR"]
C_VEN_PUL = self.parameters["circulation"]["PUL"]["C_VEN"]
C_AR_PUL = self.parameters["circulation"]["PUL"]["C_AR"]
return type(self).compute_volumes(self.parameters, self.state)

@staticmethod
def compute_volumes(parameters, state):
C_VEN_SYS = parameters["circulation"]["SYS"]["C_VEN"]
C_AR_SYS = parameters["circulation"]["SYS"]["C_AR"]
C_VEN_PUL = parameters["circulation"]["PUL"]["C_VEN"]
C_AR_PUL = parameters["circulation"]["PUL"]["C_AR"]

volumes = {
"V_LA": self.state["V_LA"],
"V_LV": self.state["V_LV"],
"V_RA": self.state["V_RA"],
"V_RV": self.state["V_RV"],
"V_AR_SYS": C_AR_SYS * self.state["p_AR_SYS"],
"V_VEN_SYS": C_VEN_SYS * self.state["p_VEN_SYS"],
"V_AR_PUL": C_AR_PUL * self.state["p_AR_PUL"],
"V_VEN_PUL": C_VEN_PUL * self.state["p_VEN_PUL"],
"V_LA": state["V_LA"],
"V_LV": state["V_LV"],
"V_RA": state["V_RA"],
"V_RV": state["V_RV"],
"V_AR_SYS": C_AR_SYS * state["p_AR_SYS"],
"V_VEN_SYS": C_VEN_SYS * state["p_VEN_SYS"],
"V_AR_PUL": C_AR_PUL * state["p_AR_PUL"],
"V_VEN_PUL": C_VEN_PUL * state["p_VEN_PUL"],
}

volumes["Heart"] = volumes["V_LA"] + volumes["V_LV"] + volumes["V_RA"] + volumes["V_RV"]
volumes["SYS"] = volumes["V_AR_SYS"] + volumes["V_VEN_SYS"]
volumes["PUL"] = volumes["V_AR_PUL"] + volumes["V_VEN_PUL"]
volumes["Total"] = volumes["Heart"] + volumes["SYS"] + volumes["PUL"]

return volumes

@property
Expand Down
16 changes: 8 additions & 8 deletions src/circulation/time_varying_elastance.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ def blanco_ventricle(
tC: float = 0.0,
TC: float = 0.4,
TR: float = 0.2,
HR: float = 1.0,
RR: float = 1.0,
):
r"""
Time-varying elastance model for the left ventricle.
Expand All @@ -24,8 +24,8 @@ def blanco_ventricle(
Duration of contraction, by default 0.4 seconds
TR : float, optional
Duration of relaxation, by default 0.2 seconds
HR : float, optional
Heart rate, by default 1.0
RR : float, optional
RR interval by default 1.0
Returns
-------
Expand Down Expand Up @@ -61,14 +61,14 @@ def blanco_ventricle(
time_rest = time_R + TR

# tC <= t <= tC + TC - Contraction
case1 = lambda t: (0 <= np.mod(t - tC, HR)) * (np.mod(t - tC, HR) < TC)
case1 = lambda t: (0 <= np.mod(t - tC, RR)) * (np.mod(t - tC, RR) < TC)
# tC + TC <= t <= tC + TC + TR - Relaxation
case2 = lambda t: (0 <= np.mod(t - time_R, HR)) * (np.mod(t - time_R, HR) < TR)
case2 = lambda t: (0 <= np.mod(t - time_R, RR)) * (np.mod(t - time_R, RR) < TR)
# tC + TC + TR <= t <= T - Rest
case3 = lambda t: 0 <= np.mod(t - time_rest, HR)
case3 = lambda t: 0 <= np.mod(t - time_rest, RR)

f_contr = lambda t: 0.5 * (1 - np.cos(np.pi / TC * (np.mod(t - tC, HR))))
f_relax = lambda t: 0.5 * (1 + np.cos(np.pi / TR * (np.mod(t - time_R, HR))))
f_contr = lambda t: 0.5 * (1 - np.cos(np.pi / TC * (np.mod(t - tC, RR))))
f_relax = lambda t: 0.5 * (1 + np.cos(np.pi / TR * (np.mod(t - time_R, RR))))
f_rest = lambda t: 0

e = lambda t: f_contr(t) * case1(t) + f_relax(t) * case2(t) + f_rest(t) * case3(t)
Expand Down
5 changes: 1 addition & 4 deletions src/circulation/zenker.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,10 @@ def default_parameters():
"R_TPR_max": 2.134 * mmHg * s / mL,
"T_sys": 4 / 15 * s,
"f_HR_min": 2 / 3 * 1 / s,
"f_HR_max": 3.0 * 1 / s,
"f_HR_max": 3 * 1 / s,
"R_valve": 0.0025 * mmHg * s / mL,
"C_PRSW_min": 25.9 * mmHg,
"C_PRSW_max": 103.8 * mmHg,
"BPM": 75.0 * units.ureg("1/minutes"),
"start_withdrawal": 0.0 * s,
"end_withdrawal": 0.0 * s,
"start_infusion": 0.0 * s,
Expand Down Expand Up @@ -166,8 +165,6 @@ def V(self, t, Pcvp, V_ES):
# REMOVE MINUS SIGN
k3 = (Pcvp + P0_LV) / R_Valve

# breakpoint()

return -(1 / k2) * np.log(
(k1 / k3) * (np.exp(-k2 * k3 * t) - 1) + np.exp(-k2 * (V_ES + k3 * t))
)
Expand Down

0 comments on commit 3b46667

Please sign in to comment.