Build your own model class

Todo

Explain how to create a custom model class

When the model is created, it must be parameterized properly. It is explained here how to create your own reference model state.

The model VanOsta2024 is build following this standard. This model is constructed by the following code.

  1"""
  2VanOsta2024 model.
  3
  4This model is first used in Van Osta et al. under review
  5"""
  6
  7from circadapt.model import Model
  8from circadapt.adapt import ModelAdapt
  9from circadapt.plot import triseg2022
 10import matplotlib.pyplot as plt
 11import numpy as np
 12
 13# define colors for the plot function
 14colors = {
 15    'Lv': np.array([213/255, 96/255, 98/255]),
 16    'Sv': np.array([236/255, 195/255, 11/255]),
 17    'Rv': np.array([6/255, 123/255, 194/255]),
 18    'Ra': np.array([27/255, 231/255, 255/255]),
 19    'La': np.array([250/255, 175/255, 150/255]),
 20    'PuArt': np.array([5/255, 75/255, 125/255]),
 21    'SyArt': np.array([150/255, 25/255, 25/255]),
 22    'RaRv': np.array([0, 0, 0]),
 23    'RvPuArt': np.array([0, 0, 0]),
 24    'LaLv': np.array([0, 0, 0]),
 25    'LvSyArt': np.array([0, 0, 0]),
 26    }
 27colors['RaRv'] = colors['Ra'] #+ 0.5*colors['Rv']
 28colors['RvPuArt'] = colors['Rv'] #+ 0.5*colors['PuArt']
 29colors['LaLv'] = colors['La'] #+ 0.5*colors['Lv']
 30colors['LvSyArt'] = colors['Lv'] #+ 0.5*colors['SyArt']
 31colors['Sv'] = (0.5*colors['Lv'] + 0.5*colors['Rv'])**1.3
 32
 33
 34class VanOsta2024(Model, ModelAdapt):
 35    def __init__(self,
 36                 solver: str = None,
 37                 path_to_circadapt: str = None,
 38                 model_state: dict = None,
 39                 ):
 40        if solver is None:
 41            solver = 'adams_moulton'
 42
 43        self._local_save_reference = True
 44
 45        ModelAdapt.__init__(self)
 46        Model.__init__(self,
 47                       solver,
 48                       path_to_circadapt=path_to_circadapt,
 49                       model_state=model_state,
 50                       )
 51
 52    def _update_settings(self):
 53        # rename model components for easy input/output
 54        self._module_rename['PressureFlowControl'] = 'PFC'
 55
 56    def build(self):
 57        pass
 58        # Circulation
 59        self.add_smart_component('ArtVen', build='SystemicCirculation')
 60        self.add_smart_component('ArtVen', build='PulmonaryCirculation')
 61        self.add_smart_component('Heart', patch_type='Patch', valve_type='Valve')
 62        self.add_smart_component('Timings')
 63        self.add_smart_component('PressureFlowControl')
 64
 65        self.set_component('PFC.circulation_volume_object', 'SyArt')
 66        self.set_component('PFC.circulation_volume_object', 'SyVen')
 67        self.set_component('PFC.circulation_volume_object', 'PuArt')
 68        self.set_component('PFC.circulation_volume_object', 'PuVen')
 69        self.set_component('PFC.circulation_volume_object', 'Peri.La')
 70        self.set_component('PFC.circulation_volume_object', 'Peri.Ra')
 71        self.set_component('PFC.circulation_volume_object', 'Peri.TriSeg.cLv')
 72        self.set_component('PFC.circulation_volume_object', 'Peri.TriSeg.cRv')
 73
 74        # # manually set papillary muscles
 75        self.set_component("Peri.RaRv.wPapMus", "Peri.TriSeg.wRv")
 76        self.set_component("Peri.LaLv.wPapMus", "Peri.TriSeg.wLv")
 77
 78    def set_reference(self):
 79        self['Chamber']['buckling'] = True
 80        self['Valve']['soft_closure'] = True
 81        self['Valve']['papillary_muscles'] = False
 82
 83        self['Solver']['dt'] = 0.001
 84        self['Solver']['dt_export'] = 0.002
 85        self.set('Solver.order',  2)
 86
 87        self.set('Model.t_cycle', 0.85)
 88        self['ArtVen']['p0'] = np.array([6306.25832487, 1000.        ])
 89        self['ArtVen']['q0'] = np.array([4.5e-05, 4.5e-05])
 90        self['ArtVen']['k'] = np.array([1., 2.])
 91        self['Tube0D']['l'] = np.array([0.4, 0.4, 0.2, 0.2])
 92        self['Tube0D']['A_wall'] = np.array([1.12362733e-04, 6.57883944e-05, 9.45910889e-05, 8.22655361e-05])
 93        self['Tube0D']['k'] = np.array([1.66666667, 2.33333333, 1.66666667, 2.33333333])
 94        self['Tube0D']['p0'] = np.array([12162.50457811,   287.7083132 ,  2132.51755623,   830.54673184])
 95        self['Tube0D']['A0'] = np.array([0.0004983 , 0.00049909, 0.00047138, 0.00050803])
 96        self['Tube0D']['target_wall_stress'] = np.array([500000., 500000., 500000., 500000.])
 97        self['Tube0D']['target_mean_flow'] = np.array([0.17, 0.17, 0.17, 0.17])
 98        self['Bag']['k'] = np.array([10.])
 99        self['Bag']['p_ref'] = np.array([1000.])
100        self['Bag']['V_ref'] = np.array([0.00054267])
101        self['Patch']['Am_ref'] = np.array([0.00425687, 0.00401573, 0.00966859, 0.00289936, 0.01084227])
102        self['Patch']['V_wall'] = np.array([4.46069398e-06, 2.14548521e-06, 7.35720515e-05, 1.88904978e-05, 3.67720116e-05,])
103        self['Patch']['v_max'] = np.array([14., 14.,  7.,  7.,  7.])
104        self['Patch']['l_se0'] = np.array([0.04, 0.04, 0.04, 0.04, 0.04])
105        self['Patch']['l_s0'] = np.array([1.8, 1.8, 1.8, 1.8, 1.8])
106        self['Patch']['l_s_ref'] = np.array([2., 2., 2., 2., 2.])
107        self['Patch']['dl_s_pas'] = np.array([0.6, 0.6, 0.6, 0.6, 0.6])
108        self['Patch']['Sf_pas'] = np.array([2248.53598   , 2684.76100348,  731.24545453,  729.06063771, 749.47694522,])
109        self['Patch']['tr'] = np.array([0.4 , 0.4 , 0.25, 0.25, 0.25])
110        self['Patch']['td'] = np.array([0.4 , 0.4 , 0.25, 0.25, 0.25])
111        self['Patch']['time_act'] = np.array([0.15 , 0.15 , 0.425, 0.425, 0.425])
112        self['Patch']['Sf_act'] = np.array([ 84000.,  84000., 120000., 120000., 120000.])
113        self['Patch']['fac_Sf_tit'] = np.array([0.01, 0.01, 0.01, 0.01, 0.01])
114        self['Patch']['k1'] = np.array([10., 10., 10., 10., 10.])
115        self['Patch']['dt'] = np.array([0., 0., 0., 0., 0.])
116        self['Patch']['C_rest'] = np.array([0., 0., 0., 0., 0.])
117        self['Patch']['l_si0'] = np.array([1.51, 1.51, 1.51, 1.51, 1.51])
118        self['Patch']['LDAD'] = np.array([1.057, 1.057, 1.057, 1.057, 1.057])
119        self['Patch']['ADO'] = np.array([0.65, 0.65, 0.65, 0.65, 0.65])
120        self['Patch']['LDCC'] = np.array([4., 4., 4., 4., 4.])
121        self['Patch']['SfPasMaxT'] = np.array([320000., 320000.,  16000.,  16000.,  16000.])
122        self['Patch']['SfPasActT'] = np.array([40000., 40000., 20000., 20000., 20000.])
123        self['Patch']['FacSfActT'] = np.array([0.8, 0.8, 1. , 1. , 1. ])
124        self['Patch']['LsPasActT'] = np.array([2.42, 2.42, 2.42, 2.42, 2.42])
125        self['Patch']['adapt_gamma'] = np.array([0.5, 0.5, 0.5, 0.5, 0.5])
126        self['Patch']['transmat00'] = np.array([-0.5751, -0.5751, -0.5751, -0.5751, -0.5751])
127        self['Patch']['transmat01'] = np.array([-0.7851, -0.7851, -0.7851, -0.7851, -0.7851])
128        self['Patch']['transmat02'] = np.array([0.6063, 0.6063, 0.6063, 0.6063, 0.6063])
129        self['Patch']['transmat03'] = np.array([-0.5565, -0.5565, -0.5565, -0.5565, -0.5565])
130        self['Patch']['transmat10'] = np.array([-0.1279, -0.1279, -0.1279, -0.1279, -0.1279])
131        self['Patch']['transmat11'] = np.array([0.0999, 0.0999, 0.0999, 0.0999, 0.0999])
132        self['Patch']['transmat12'] = np.array([0.2066, 0.2066, 0.2066, 0.2066, 0.2066])
133        self['Patch']['transmat13'] = np.array([-1.8441, -1.8441, -1.8441, -1.8441, -1.8441])
134        self['Patch']['transmat20'] = np.array([-0.1865, -0.1865, -0.1865, -0.1865, -0.1865])
135        self['Patch']['transmat21'] = np.array([-0.201, -0.201, -0.201, -0.201, -0.201])
136        self['Patch']['transmat22'] = np.array([1.3195, 1.3195, 1.3195, 1.3195, 1.3195])
137        self['Patch']['transmat23'] = np.array([-11.8745, -11.8745, -11.8745, -11.8745, -11.8745])
138        self['Valve']['adaptation_A_open_fac'] = np.array([1., 1., 1., 1., 1., 1.])
139        self['Valve']['A_open'] = np.array([0.00049916, 0.00047155, 0.00047155, 0.00050805, 0.00049835,
140               0.00049835])
141        self['Valve']['A_leak'] = np.array([0.00049916, 1.e-09, 1.e-09, 0.00050805, 1.e-09, 1.e-09])
142        self['Valve']['l'] = np.array([0.01260512, 0.01225144, 0.01225144, 0.01271679, 0.01259479,
143               0.01259479])
144        self['Valve']['L_fac_prox'] = np.array([0.75, 0.75, 0.75, 0.75, 0.75, 0.75])
145        self['Valve']['L_fac_dist'] = np.array([0.75, 0.75, 0.75, 0.75, 0.75, 0.75])
146        self['Valve']['L_fac_valve'] = np.array([1.5, 1.5, 1.5, 1.5, 1.5, 1.5])
147        self['Valve']['rho_b'] = np.array([1050., 1050., 1050., 1050., 1050., 1050.])
148        self['Valve']['papillary_muscles'] = np.array([False, False, False, False, False, False])
149        self['Valve']['papillary_muscles_slope'] = np.array([100., 100., 100., 100., 100., 100.])
150        self['Valve']['papillary_muscles_min'] = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
151        self['Valve']['papillary_muscles_A_open_fac'] = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
152        self['Valve']['soft_closure'] = np.array([ True,  True,  True,  True,  True,  True])
153        self['Valve']['fraction_A_open_Aext'] = np.array([0.9, 0.9, 0.9, 0.9, 0.9, 0.9])
154        self['Timings']['time_fac'] = np.array([1.])
155        self['Timings']['tau_av'] = np.array([0.150025])
156        self['Timings']['dtau_av'] = np.array([0.])
157        self['Timings']['law_tau_av'] = np.array([1])
158        self['Timings']['law_Ra2La'] = np.array([1])
159        self['Timings']['law_ta'] = np.array([1])
160        self['Timings']['law_tv'] = np.array([1])
161        self['Timings']['c_tau_av0'] = np.array([0.])
162        self['Timings']['c_tau_av1'] = np.array([0.1765])
163        self['Timings']['c_ta_rest'] = np.array([0.])
164        self['Timings']['c_ta_tcycle'] = np.array([0.17647059])
165        self['Timings']['c_tv_rest'] = np.array([0.1])
166        self['Timings']['c_tv_tcycle'] = np.array([0.4])
167        self['PFC']['p0'] = np.array([12200.])
168        self['PFC']['q0'] = np.array([8.5e-05])
169        self['PFC']['stable_threshold'] = np.array([0.0001])
170        self['PFC']['is_active'] = np.array([ True])
171        self['PFC']['fac'] = np.array([1.])
172        self['PFC']['alpha'] = np.array([1.])
173        self['PFC']['epsilon'] = np.array([0.4])
174        self['TriSeg']['tau'] = np.array([2.])
175        self['TriSeg']['ratio_septal_LV_Am'] = np.array([0.3])
176        self['TriSeg']['max_number_of_iterations'] = np.array([100.])
177        self['TriSeg']['thresh_F'] = np.array([0.001])
178        self['TriSeg']['thresh_dV'] = np.array([1.e-09])
179        self['TriSeg']['thresh_dY'] = np.array([1.e-06])
180        self['TriSeg']['collapse_eps'] = 0.0
181        self['TriSeg']['collapse_start'] = 0.25
182
183        self.run(stable=True)
184        # self.adapt()
185
186        # set target volume
187        self['PFC']['target_volume'] = self['PFC']['circulation_volume']
188
189        return
190
191
192    def get_unittest_targets(self):
193        """Hardcoded results after initializing and running 1 beat."""
194        return {
195            'LVEDV': 120.3,
196            'LVESV':  48.0,
197            }
198
199    def get_unittest_results(self, model):
200        """Real-time results after initializing and running 1 beat."""
201        LVEDV = np.max(model['Cavity']['V'][:, 'cLv'])*1e6
202        LVESV = np.min(model['Cavity']['V'][:, 'cLv'])*1e6
203        return {
204            'LVEDV': LVEDV,
205            'LVESV': LVESV,
206            }
207
208
209    def plot(self, fig=None):
210        # TODO
211        self.plot_extended(fig)
212
213    def plot_extended(self, fig=None):
214        if len(self['Solver']['t'])==0:
215            raise RuntimeError('Can not plot the model, as no data is available. Did you simulate the model?')
216
217        if fig is None:
218            fig = 1
219        if isinstance(fig, int):
220            fig = plt.figure(fig, clear=True, figsize=(12, 8))
221
222        # Settings
223        grid_size = [32, 32]
224
225        def get_lim(module, signal, locs=slice(None, None, None)):
226            signal = self[module][signal][:, locs]
227            lim = np.array([np.min(signal), np.max(signal)])
228            lim += np.array([-1, 1]) * 0.1*np.diff(lim)
229            return lim
230
231        lim_V = get_lim('Cavity', 'V', ['cRv', 'Ra', 'La', 'cLv']) * 1e6
232        lim_V[0] = 0
233        lim_p = get_lim('Cavity', 'p', ['cRv', 'Ra', 'La', 'cLv']) / 133
234        lim_p[0] = np.min([lim_p[0], 0])
235        lim_Ls = get_lim('Patch', 'l_s')
236        lim_Sf = get_lim('Patch', 'Sf') * 1e-3
237        lim_q = get_lim('Valve', 'q',
238                        ['LaLv', 'RaRv', 'LvSyArt', 'RaPuArt']) * 1e6
239
240        all_lim = [lim_V, lim_p, lim_Ls, lim_Sf, lim_q]
241        if (np.any(np.isnan(all_lim)) or np.any(np.isinf(all_lim))):
242            lim_V = [0, 200]
243            lim_p = [0, 150]
244            lim_Ls = [1.5, 2.0]
245            lim_Sf = [0, 100]
246            lim_q = [-1e-3, 1e-3]
247
248        # Pressure Volume plot
249        axPV = plt.subplot2grid(grid_size, (0, 17), rowspan=15, colspan=15, fig=fig)
250        axPV.plot(self['Cavity']['V'][:, 'cLv']*1e6,
251                  self['Cavity']['p'][:, 'cLv']/133,
252                  c=colors['Lv'],
253                  )
254        axPV.plot(self['Cavity']['V'][:, 'cRv']*1e6,
255                  self['Cavity']['p'][:, 'cRv']/133,
256                  c=colors['Rv'],
257                  )
258        axPV.plot(self['Cavity']['V'][:, 'La']*1e6,
259                  self['Cavity']['p'][:, 'La']/133,
260                  c=colors['La'],
261                  )
262        axPV.plot(self['Cavity']['V'][:, 'Ra']*1e6,
263                  self['Cavity']['p'][:, 'Ra']/133,
264                  c=colors['Ra'],
265                  )
266        axPV.spines[['top', 'right']].set_visible(False)
267        axPV.set_title('Pressure-Volume loop', weight='bold')
268        axPV.set_xlabel('Volume [mL]')
269        axPV.set_ylabel('Pressure [mmHg]')
270        axPV.spines[['bottom', 'left']].set_position(('outward', 5))
271
272        ylabel_x_left = -0.25
273        ylabel_x_right = 1.25
274
275        # Volumes
276        t = self['Solver']['t']*1e3
277
278        axVRv = plt.subplot2grid(grid_size, (0, 0), rowspan=8, colspan=6, fig=fig)
279        axVRv.plot(t, self['Cavity']['V'][:, 'cRv']*1e6,
280                   c=colors['Rv'],
281                   )
282        axVRv.plot(t, self['Cavity']['V'][:, 'Ra']*1e6,
283                   c=colors['Ra'],
284                   )
285        axVRv.set_ylim(lim_V)
286        axVRv.set_ylabel('Volume\n[mL]')
287        axVRv.spines[['top', 'right']].set_visible(False)
288        axVRv.set_title('Right Heart', weight='bold')
289        # axVRv.set_xticks([])
290        axVRv.tick_params(axis='both', direction='in')
291        axVRv.yaxis.set_label_coords(ylabel_x_left, 0.5)
292
293
294        axVLv = plt.subplot2grid(grid_size, (0, 6), rowspan=8, colspan=6, fig=fig)
295        axVLv.plot(t, self['Cavity']['V'][:, 'cLv']*1e6,
296                   c=colors['Lv'],
297                   )
298        axVLv.plot(t, self['Cavity']['V'][:, 'La']*1e6,
299                   c=colors['La'],
300                   )
301        axVLv.set_ylabel('Volume\n[mL]')
302        axVLv.set_ylim(lim_V)
303        axVLv.yaxis.set_ticks_position('right')
304        axVLv.yaxis.set_label_position('right')
305        axVLv.spines['right'].set_position(('outward', 0))
306        axVLv.spines[['top', 'left']].set_visible(False)
307        axVLv.set_title('Left Heart', weight='bold')
308        # axVLv.set_xticks([])
309        axVLv.tick_params(axis='both', direction='in')
310        axVLv.yaxis.set_label_coords(ylabel_x_right, 0.5)
311
312        # Pressures
313        axpRv = plt.subplot2grid(
314            grid_size, (8, 0), rowspan=8, colspan=6, fig=fig, sharex=axVRv)
315        axpRv.plot(t, self['Cavity']['p'][:, 'cRv']/133,
316                   c=colors['Rv'],
317                   )
318        axpRv.plot(t, self['Cavity']['p'][:, 'Ra']/133,
319                   c=colors['Ra'],
320                   )
321        axpRv.plot(t, self['Cavity']['p'][:, 'PuArt']/133,
322                   c=colors['PuArt'],
323                   )
324        axpRv.spines[['top', 'right']].set_visible(False)
325        # axpRv.set_xticks([])
326        axpRv.tick_params(axis='both', direction='in')
327        axpRv.set_ylim(lim_p)
328        axpRv.set_ylabel('Pressure\n[mmHg]')
329        axpRv.yaxis.set_label_coords(ylabel_x_left, 0.5)
330
331        axpLv = plt.subplot2grid(grid_size, (8, 6), rowspan=8, colspan=6, fig=fig, sharex=axVLv)
332        axpLv.plot(t, self['Cavity']['p'][:, 'cLv']/133,
333                   c=colors['Lv'],
334                   )
335        axpLv.plot(t, self['Cavity']['p'][:, 'La']/133,
336                   c=colors['La'],
337                   )
338        axpLv.plot(t, self['Cavity']['p'][:, 'SyArt']/133,
339                   c=colors['SyArt'],
340                   )
341        axpLv.yaxis.set_ticks_position('right')
342        axpLv.yaxis.set_label_position('right')
343        axpLv.spines['right'].set_position(('outward', 0))
344        axpLv.spines[['top', 'left']].set_visible(False)
345        # axpLv.set_xticks([])
346        axpLv.tick_params(axis='both', direction='in')
347        axpLv.set_ylim(lim_p)
348        axpLv.set_ylabel('Pressure\n[mmHg]')
349        axpLv.yaxis.set_label_coords(ylabel_x_right, 0.5)
350
351        # Valves
352        ax = plt.subplot2grid(grid_size, (16, 0), rowspan=6, colspan=6, fig=fig, sharex=axVRv)
353        ax.plot(t, self['Valve']['q'][:, 'RaRv']*1e6,
354                c=colors['RaRv'],
355                )
356        ax.plot(t, self['Valve']['q'][:, 'RvPuArt']*1e6,
357                c=colors['RvPuArt'],
358                )
359        ax.spines[['top', 'right']].set_visible(False)
360        ax.set_ylim(lim_q)
361        # ax.set_xticks([])
362        ax.set_ylabel('Flow\n[mL/s]')
363        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
364
365        ax = plt.subplot2grid(grid_size, (16, 6), rowspan=6, colspan=6, fig=fig, sharex=axVLv)
366        ax.plot(t, self['Valve']['q'][:, 'LaLv']*1e6,
367                   c=colors['LaLv'],
368                   )
369        ax.plot(t, self['Valve']['q'][:, 'LvSyArt']*1e6,
370                   c=colors['LvSyArt'],
371                   )
372        ax.spines[['top', 'left']].set_visible(False)
373        ax.set_ylim(lim_q)
374        ax.yaxis.set_ticks_position('right')
375        ax.yaxis.set_label_position('right')
376        # ax.set_xticks([])
377        ax.set_ylabel('Flow\n[mL/s]')
378        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
379
380        # Stress
381        ax = plt.subplot2grid(grid_size, (22, 0), rowspan=4, colspan=6, fig=fig, sharex=axVRv)
382        ax.plot(t, self['Patch']['Sf'][:, 'pRv0']*1e-3,
383                   c=colors['Rv'],
384                   )
385        ax.plot(t, self['Patch']['Sf'][:, 'pRa0']*1e-3,
386                   c=colors['Ra'],
387                   )
388        ax.spines[['top', 'right']].set_visible(False)
389        # ax.set_xticks([])
390        ax.set_ylim(lim_Sf)
391        ax.set_ylabel('Total\nstress [kPa]')
392        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
393
394        ax = plt.subplot2grid(grid_size, (22, 6), rowspan=4, colspan=6, fig=fig, sharex=axVLv)
395        ax.plot(t, self['Patch']['Sf'][:, 'pLv0']*1e-3,
396                   c=colors['Lv'],
397                   )
398        ax.plot(t, self['Patch']['Sf'][:, 'pSv0']*1e-3,
399                   c=colors['Sv'],
400                   )
401        ax.plot(t, self['Patch']['Sf'][:, 'pLa0']*1e-3,
402                   c=colors['La'],
403                   )
404        ax.spines[['top', 'left']].set_visible(False)
405        ax.yaxis.set_ticks_position('right')
406        ax.yaxis.set_label_position('right')
407        # ax.set_xticks([])
408        ax.set_ylim(lim_Sf)
409        ax.set_ylabel('Total\nstress [kPa]')
410        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
411
412        # Sarcomere Length
413        ax = plt.subplot2grid(grid_size, (26, 0), rowspan=6, colspan=6, fig=fig, sharex=axVRv)
414        ax.plot(t, self['Patch']['l_s'][:, 'pRv0'],
415                   c=colors['Rv'],
416                   )
417        ax.plot(t, self['Patch']['l_s'][:, 'pRa0'],
418                   c=colors['Ra'],
419                   )
420        ax.spines[['top', 'right']].set_visible(False)
421        ax.spines['bottom'].set_position(('outward', 5))
422        ax.set_ylim(lim_Ls)
423        ax.set_ylabel('Sarcomere\nlength [$\\mu$m]')
424        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
425
426        ax = plt.subplot2grid(grid_size, (26, 6), rowspan=6, colspan=6, fig=fig, sharex=axVLv)
427        ax.plot(t, self['Patch']['l_s'][:, 'pLv0'],
428                   c=colors['Lv'],
429                   )
430        ax.plot(t, self['Patch']['l_s'][:, 'pSv0'],
431                   c=colors['Sv'],
432                   )
433        ax.plot(t, self['Patch']['l_s'][:, 'pLa0'],
434                   c=colors['La'],
435                   )
436        ax.spines[['top', 'left']].set_visible(False)
437        ax.spines['bottom'].set_position(('outward', 5))
438        ax.yaxis.set_ticks_position('right')
439        ax.yaxis.set_label_position('right')
440        ax.set_ylim(lim_Ls)
441        ax.set_ylabel('Sarcomere\nlength [$\\mu$m]')
442        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
443
444        # ax.set_xlabel('Time [ms]')
445        # ax.xaxis.set_label_coords(0, -0.3)
446
447
448        # Plot TriSeg
449        titles = ['Pre-A', 'Onset QRS', 'Peak LV \n pressure', 'AV close']
450        idx = [0,
451               np.argmax(np.diff(self['Patch']['C'][:, 'pLv0'])>0),
452               np.argmax(self['Cavity']['p'][:, 'cLv']),
453               len(t) - 1 - np.argmax(
454                   np.diff(self['Valve']['q'][:, 'LvSyArt'][::-1])>0)
455               ]
456
457        for i in range(4):
458            ax = plt.subplot2grid(grid_size, (26, 16+4*i),
459                                  rowspan=5, colspan=4, fig=fig)
460            triseg2022(self, ax, idx[i],
461                       colors=[colors['Lv'], colors['Sv'], colors['Rv']])
462            ax.spines[['top', 'right', 'bottom', 'left']].set_visible(False)
463            ax.set_xticks([])
464            ax.set_yticks([])
465            plt.xlabel(titles[i], fontsize=12)
466
467        # Plot settings
468        plt.subplots_adjust(
469            top=0.96,
470            bottom=0.05,
471            left=0.075,
472            right=0.98,
473            hspace=5,
474            wspace=0.5)
475        plt.draw()
476
477        # Stress strain
478        ax_right_stress_strain = plt.subplot2grid(grid_size, (18, 17),
479                              rowspan=7, colspan=7, fig=fig)
480        ax_left_stress_strain = plt.subplot2grid(grid_size, (18, 24),
481                              rowspan=7, colspan=7, fig=fig)
482
483
484        ax_right_stress_strain.plot(
485            self['Patch']['Ef'][:, 'pRa0'],
486            self['Patch']['Sf'][:, 'pRa0']*1e-3,
487            c=colors['Ra'],
488            )
489        ax_right_stress_strain.plot(
490            self['Patch']['Ef'][:, 'pRv0'],
491            self['Patch']['Sf'][:, 'pRv0']*1e-3,
492            c=colors['Rv'],
493            )
494        ax_left_stress_strain.plot(
495            self['Patch']['Ef'][:, 'pLa0'],
496            self['Patch']['Sf'][:, 'pLa0']*1e-3,
497            c=colors['La'],
498            )
499        ax_left_stress_strain.plot(
500            self['Patch']['Ef'][:, 'pLv0'],
501            self['Patch']['Sf'][:, 'pLv0']*1e-3,
502            c=colors['Lv'],
503            )
504        ax_left_stress_strain.plot(
505            self['Patch']['Ef'][:, 'pSv0'],
506            self['Patch']['Sf'][:, 'pSv0']*1e-3,
507            c=colors['Sv'],
508            )
509
510        ylim = [np.min(self['Patch']['Sf'])*1e-3,
511                np.max(self['Patch']['Sf'])*1e-3]
512        ylim_ptp = np.ptp(ylim)*0.025
513        ylim[0] -= ylim_ptp
514        ylim[1] += ylim_ptp
515        ax_right_stress_strain.set_ylim(ylim)
516        ax_left_stress_strain.set_ylim(ylim)
517        xlim = [np.min(self['Patch']['Ef']),
518                np.max(self['Patch']['Ef'])]
519        xlim_ptp = np.ptp(xlim)*0.025
520        xlim[0] -= xlim_ptp
521        xlim[1] += xlim_ptp
522        ax_right_stress_strain.set_xlim(xlim)
523        ax_left_stress_strain.set_xlim(xlim)
524
525        ax_right_stress_strain.spines[['left', 'bottom']].set_position(('outward', 5))
526        ax_left_stress_strain.spines[['right', 'bottom']].set_position(('outward', 5))
527
528        ax_right_stress_strain.spines[['top', 'right']].set_visible(False)
529        ax_right_stress_strain.set_ylabel('Stress\n[kPa]', rotation=0, va='bottom', ha='right')
530        ax_right_stress_strain.yaxis.set_label_coords(0., 1.05)
531        ax_right_stress_strain.set_xlabel('Strain [-]')
532        ax_left_stress_strain.spines[['top', 'left']].set_visible(False)
533        ax_left_stress_strain.yaxis.set_ticks_position('right')
534        ax_left_stress_strain.yaxis.set_label_position('right')
535        ax_left_stress_strain.set_ylabel('Stress\n[kPa]', rotation=0, va='bottom', ha='left')
536        ax_left_stress_strain.yaxis.set_label_coords(1., 1.05)
537        ax_left_stress_strain.set_xlabel('Strain [-]')
538