Preload Afterload experiment

Tags: PreAfterloadExperiment Patch2022

Tutorial on how to implement the preload-afterload experiment. This tutorial discusses how to build the model, how to switch between different patch types, and how to extract and plot information from the model.

Todo

Make Preload Afterload methods figure

As all projects, we start with importing packages.

 1# sys.path.append('../../../src/')
 2
 3import circadapt
 4# Uncomment next lines if not installed
 5# circadapt.DEFAULT_PATH_TO_CIRCADAPT = "../../../CircAdapt_Library/out/build/x64-Release/CircAdaptLib.dll"
 6
 7
 8
 9# import
10from circadapt import CircAdapt
11import numpy as np
12import matplotlib.pyplot as plt
13
14import time
15

We then create an empty model. Then, we add the PreAfterloadExperiment wall object. This object is empty by default, but it needs at least 1 patch object to work. In this tutorial, we use the Patch2022 object.

1model = CircAdapt('forward_euler')
2model.add_component('PreAfterloadExperiment', 'PAE')
3model.add_component('Patch2022', 'P', 'PAE')
4

The model has to be parametized, as default parameter values do not make sense for this setup.

 1model.set('Model.t_cycle', 1e-0)
 2model.set('Solver.dt', 1e-3)
 3model.set('Solver.dt_export', 1e-3)
 4
 5model['PreAfterloadExperiment']['Am_ref_afterload'] =  0.006
 6model['PreAfterloadExperiment']['T_afterload'] =  200
 7model['PreAfterloadExperiment']['n_iter'] =  5
 8model['Patch2022']['Am_ref'] =  0.005
 9model['Patch2022']['V_wall'] = 92.43*1e-6
10model['Patch2022']['ADO'] = 3.0
11model['Patch2022']['tr'] = 0.5
12model['Patch2022']['td'] = 0.5
13model['Patch2022']['k1'] = 10.
14model['Patch2022']['v_max'] = 7.0
15model['Patch2022']['dt'] = 0.1
16
17model.set('Model.PAE.P.l_si', -0.04 + 2 * np.sqrt(
18    model['PreAfterloadExperiment']['Am_ref_afterload'][0]/model['Patch2022']['Am_ref'][0]))
19

Now the model setup is finished and we can run the simulation. For this setup, only 1 run is sufficient.

1t0 = time.time()
2model.run(1)
3t1 = time.time()
4print(t1-t0)
5

The simulation will result in the stresses and tensions shown in the plot below.

(Source code, png, hires.png, pdf)

../../../_images/Patch2022.png

The full code is shown below.

 1# sys.path.append('../../../src/')
 2
 3import circadapt
 4# Uncomment next lines if not installed
 5# circadapt.DEFAULT_PATH_TO_CIRCADAPT = "../../../CircAdapt_Library/out/build/x64-Release/CircAdaptLib.dll"
 6
 7
 8
 9# import
10from circadapt import CircAdapt
11import numpy as np
12import matplotlib.pyplot as plt
13
14import time
15
16# %% Create custom model
17model = CircAdapt('forward_euler')
18model.add_component('PreAfterloadExperiment', 'PAE')
19model.add_component('Patch2022', 'P', 'PAE')
20
21# %% Set model parameters
22model.set('Model.t_cycle', 1e-0)
23model.set('Solver.dt', 1e-3)
24model.set('Solver.dt_export', 1e-3)
25
26model['PreAfterloadExperiment']['Am_ref_afterload'] =  0.006
27model['PreAfterloadExperiment']['T_afterload'] =  200
28model['PreAfterloadExperiment']['n_iter'] =  5
29model['Patch2022']['Am_ref'] =  0.005
30model['Patch2022']['V_wall'] = 92.43*1e-6
31model['Patch2022']['ADO'] = 3.0
32model['Patch2022']['tr'] = 0.5
33model['Patch2022']['td'] = 0.5
34model['Patch2022']['k1'] = 10.
35model['Patch2022']['v_max'] = 7.0
36model['Patch2022']['dt'] = 0.1
37
38model.set('Model.PAE.P.l_si', -0.04 + 2 * np.sqrt(
39    model['PreAfterloadExperiment']['Am_ref_afterload'][0]/model['Patch2022']['Am_ref'][0]))
40
41# %% Run model
42t0 = time.time()
43model.run(1)
44t1 = time.time()
45print(t1-t0)
46
47# %% Plot model
48fig = plt.figure(1, clear=True, figsize=(12, 8))
49m = 2
50n = 3
51
52t = model.get('Solver.t') * 1e3
53
54ax1 = fig.add_subplot(m, n, 1)
55ax1.plot(t, model.get('Model.PAE.P.l_s'))
56ax1.plot(t, model.get('Model.PAE.P.l_si'))
57ax1.legend(['l_s', 'l_si'])
58ax1.set_ylabel('Sarcomere length [$\mu m$]')
59
60ax4 = fig.add_subplot(m, n, 4)
61ax4.plot(t, model.get('Model.PAE.P.Am')*1e4, label='Patch Am')
62ax4.plot(t, model.get('Model.PAE.P.Am0')*1e4, label='Patch Am0')
63ax4.set_ylabel('Segment area [$cm^2$]')
64ax4.legend()
65
66ax3 = fig.add_subplot(m, n, 3)
67ax3.plot(t, model.get('Model.PAE.P.Sf')*1e-3, label='Sf')
68ax3.plot(t, model.get('Model.PAE.P.Sf_ecm')*1e-3, label='SfEcm')
69ax3.set_ylabel('Stress [kPa]')
70ax3.legend()
71
72ax2 = fig.add_subplot(m, n, 2)
73ax2.plot(t, model.get('Model.PAE.P.C'), label='C')
74ax2.set_ylabel('C-curve [-]')
75
76ax5 = fig.add_subplot(m, n, 5)
77ax5.plot(t, model.get('Model.PAE.dA_dT')*1e4, label='Tension Wall')
78ax5.set_ylabel('dA_dT [$cm^2/N$]')
79
80ax6 = fig.add_subplot(m, n, 6)
81ax6.plot(t, model.get('Model.PAE.T'), label='Tension Wall')
82ax6.set_ylabel('Tension [N]')
83
84# Plot design
85for ax in [ax1, ax2, ax3, ax4, ax5, ax6]:
86    ax.spines[['right', 'top']].set_visible(False)
87ax5.set_xlabel('Time [ms]')
88
89plt.tight_layout()