87 The origin of third heart sound - II

/*----Divider----*/
The model in this episode is based on this article:
Annals of Biomedical Engineering, Vol. 19, pp. 651-667, 1991.
/*----Divider----*/

In our last episode we have introduced a schematic heart in diastolic phase, and from basic pressure-flow relationship, the pressure gradient between pulmonary veins and left ventricle during diastolic phase could be stated as below.

$P_{PV}=P_{LV}+Q_{MV}R_{PV}$

Since $Q_{MV}$ is related to the volume change of left ventricle, if we could properly relate $V_{LV}$ to $P_{LV}$, the differential equation could be solved. First we will look at the elastic part of $P_{LV}$. Out of instinct, we could express it as follow:

$P_{c} = b(V_{LV}-V_{0})$

This is based on the assumption of a constant stress-stiffness relationship as stated by Hooke's law. However, experiments on isolated heart muscle showed that the stress-stiffness is in fact linear[1]:

$\frac{d\sigma}{d\epsilon} = k\sigma+c$

This implies an exponential stress-strain response as you solve the differential equation:

$\sigma = \frac{c}{k}(e^{k\epsilon}-1)$

Therefore, in our model, we will assume the elastic pressure of left ventricle to be:

$P_{c} = b_{1}(e^{b_{2}(V_{LV} - V_{0})} - 1)$

There is of course an active component of ventricular pressure when it is in systolic phase and early diastolic phase:

$P_{o} = \frac{1}{2}K( 1 - \cos(\frac{2\pi t}{T_{s}}))(V_{LV} - V_{d})$

Note that generally $V_{d}$ is not equal to $V_{0}$ since the former stands for the volume at which no active pressure could be generated by ventricle, or the dead space volume, while the latter means the volume where there is no passive wall stress.

The last part of our model is the viscous response of ventricular wall, $P_{v}$ (i.e., the ventricular pressure generated by how fast its volume changes), which can be formulated as Maxwell model:

$\frac{dV_{LV}}{dt} = \frac{P_{v}}{R_{LV}} + C_{LV}\frac{dP_{v}}{dt}$

The whole model of our heart could be summarized as the Figure 1 of the original article. The ventricular pressure is composed of elastic, viscous and active components. Pulmonary vein is a constant pressure source for venous return, and aorta is a compliant reservoir which provides incessant blood flow in both systolic and diastolic phase.

Figure 1 of the original article (Annals of Biomedical Engineering, Vol. 19, pp. 651-667, 1991.) The figure summarized different components of the heart model. The ventricular pressure is composed of 3 components -- elastic pressure, viscous pressure and active pressure. The venous return from pulmonary veins is a constant pressure source, and aorta acts as a reservoir to provide incessant blood flow. Mitral valves and aortic valves act as a diode to prevent retrograde blood flow.

From this model, we could simulate the ventricular pressure in a cardiac cycle. The heart sound is determined by filtering the pressure to obtain pressure waves between 25 and 160 Hz. With all these prepared, we are ready to have our Python code.

from scipy.integrate import odeint
import numpy as np
from numpy import arange
from pylab import *
from scipy.signal import iirfilter, filtfilt, butter
import matplotlib.pyplot as plt

class Heart:
 def __init__(self, HR=0, Ppv=0, Cs=0, Ro=0, b1=0, b2=0, Vo=0, Rlv=0, Clv=0, K=0, Vd=0):
  self.HR=HR
  self.Tcycle = 60/HR
  self.Ts = 0.32
  self.Ppv = Ppv
  self.Rpv = 0.04
  self.Rs = 0.90
  self.Cs = Cs
  self.Ro = Ro
  self.b1 = b1
  self.b2 = b2
  self.Vo = Vo
  self.Rlv = Rlv
  self.Clv = Clv
  self.K = K
  self.Vd = Vd
def heartModel(state, t, H, dt):
 Vlv, Pv, Pao = state
 Po = 1/2*H.K*(1-np.cos(2*pi*np.mod(t, H.Tcycle)/H.Ts))*(Vlv - H.Vd)*(np.mod(t, H.Tcycle) <= H.Ts)
 Pc = H.b1*(np.exp(H.b2*(Vlv - H.Vo)) - 1)
 Plv = Po + Pc + Pv
 Qmv = (H.Ppv - Plv)/H.Rpv*(H.Ppv > Plv)
 Qao = (Plv - Pao)/H.Ro*(Plv > Pao)
 d_Pao = (Qao - Pao/H.Rs)/H.Cs
 d_Vlv = (Qmv - Qao)
 d_Pv =  (Qmv - Qao - Pv/H.Rlv)/H.Clv
 return [d_Vlv, d_Pv, d_Pao]
def calHeart(state, H):
    Vlv = state[:, 0]
    Pv = state[:, 1]
    Po = 1/2*H.K*(1-np.cos(2*pi*np.mod(t, H.Tcycle)/H.Ts))*(Vlv - H.Vd)*(np.mod(t, H.Tcycle) <= H.Ts)
    Pc = H.b1*(np.exp(H.b2*(Vlv - H.Vo)) - 1)
    Plv = Po + Pc + Pv
    return Po, Pc, Plv
def fourierHeart(P):
    b, a = iirfilter(4, [0.05, 0.32], 1, 60, analog=False, ftype='cheby1')
    filtState = filtfilt(b, a, P)
    b, a = butter(3, 0.20, analog=True)
    filtState2 = filtfilt(b, a, filtState)
    return filtState2
def plotHeart(state, H, TEXT):
    Po, Pc, Plv = calHeart(state, H)
    Pao = state[:,2]
    Pv = state[:,1]
    Vlv = state[:,0]
    filtState = fourierHeart(Plv)
    #pressure changes curve
    plt.figure()
    plt.plot(t, Plv, 'b-', t, Pv, 'm-', t, Pc, 'y-', t, Po, 'g-', t, Pao, 'r-')
    plt.xlabel('time (sec)', fontsize=18)
    plt.ylabel('pressure (mmHg)', fontsize=18) 
    plt.legend(('LV pressure', 'viscous pressure', 'elastic pressure', 'active pressure', 'Ao pressure'), loc = 'best', fontsize=16)
    plt.title('pressure changes during a cardiac cycle, ' +TEXT, fontsize=18)
    plt.grid(True)
    plt.savefig('Pressure_'+TEXT+'.png', dpi=150)
    plt.show()
    #pressure & sound
    plt.figure()
    fig, ax1 = plt.subplots()
    ax1.plot(t, Plv, 'b-',  t, Pv, 'm-', t, Pc, 'y-', t, Po, 'g-')
    ax1.set_xlabel('time (s)', fontsize=18)
    ax1.set_ylabel('LV pressure(mmHg)', color='b', fontsize=18)
    ax1.tick_params('y', colors='b')
    ax2 = ax1.twinx()
    ax2.plot(t, filtState, 'r')    
    ax2.set_ylabel('heart sound', color='r', fontsize=18)
    ax2.tick_params('y', colors='r')
    fig.tight_layout()
    title('heart sound, ' + TEXT, fontsize=14)
    ax2.plot(np.nan, 'b-', np.nan, 'm-', np.nan, 'y-', np.nan, 'g-')
    ax2.legend(('heart sound', 'LV pressure', 'viscous pressure', 'elastic pressure', 'active pressure'), loc = 'best', fontsize=16)
    savefig('heart sound_'+TEXT+'.png')
    show()
    #pressure & volume
    figure()
    fig, ax1 = plt.subplots()
    ax1.plot(t, Plv, 'b-',  t, Pv, 'm-', t, Pc, 'y-', t, Po, 'g-')
    ax1.set_xlabel('time (s)', fontsize=18)
    ax1.set_ylabel('LV pressure(mmHg)', color='b', fontsize=18)
    ax1.tick_params('y', colors='b')
    ax2 = ax1.twinx()
    ax2.plot(t, Vlv, 'r')
    ax2.set_ylabel('LV volume(mL)', color='r', fontsize=18)
    ax2.plot(np.nan, 'b-', np.nan, 'm-', np.nan, 'y-', np.nan, 'g-')
    ax2.legend(('LV volume', 'LV pressure', 'viscous pressure', 'elastic pressure', 'active pressure'), loc = 'best', fontsize=16)
    ax2.tick_params('y', colors='r')
    fig.tight_layout()
    title('LV pressure & volume changes, ' + TEXT, fontsize=14)
    savefig('PV_change_'+TEXT+'.png')
    plt.show()

dt = 0.001
t = arange(0, 2, dt)

H_normal = Heart(80.0, 8.0, 7.00, 0.02, 5.00, 0.008, 40.00, 0.04, 2.00, 5.00, 2.00)
init_state_normal = [85, 0, 80]
state_normal = odeint(heartModel, init_state_normal, t, args=(H_normal, dt))
plotHeart(state_normal, H_normal, 'normal_heart')

H_DCM = Heart(90.0, 25.0, 3.00, 0.05, 3.00, 0.004, 10.00, 5.00, 0.30, 1.00, -10.00)
init_state_DCM = [250, 0, 80]
state_DCM = odeint(heartModel, init_state_DCM, t, args = (H_DCM, dt))
plotHeart(state_DCM, H_DCM, 'DCM') 

And the code will generate some figures like these:

We will examine these figures more carefully in our next episode.

/*----Reference----*/
[1] Circ Res. 1973 Aug; 33(2):233-43.