Question

3. Write a Matlab script that describes the dynamics of Pressure and Flow in the systemic arteries and Left Ventricle. Assume that :

1) Pressure in the Left atrium (P_LA) and in the Systemic Veins (Psv) remain constant.

2) A time dependent LV compliance

3) Opening/Closing of the heart valves instantaneously with the direction of flow (i.e. valves are at an open or closed states. Use parameter values from the handout.

Answer 1

/* define normal human, hemodynamic constants in 2-element Windkessel model */ pl = 72    /* cardiac cycles per minute */ s = 60/pl       /* period of cardiac cycle in seconds */ h = (2/5)*s    /* period of systole in seconds */ ml = 90      /* blood output to aorta or pulmonary artery per cardiac cycle in cm^3 */ /* define function for blood flow as variable amplitude sine wave */ fct i(t,q0) = q0*sin(pi*t/h) /* find peak amplitude of flow such that the total cardiac output is ml */ fct ampl() = root(z,0,1000,ml-integral(t,0,h,i(t,z))) q0 = ampl() /* define and compute flow for 3 cardiac cycles */ fct i1(t) = if ((tmd_mod(t,s)) 

The flow function, aortic pressure, and pulmonary pressure curves can be graphed with the following additional MLAB commands:

draw af left title "Blood Flow (ml/sec)" size .007 frame 0 to 1, .66 to .99 w1 = w draw ms left title "Aortic Pressure (mmHg)" size .007 frame 0 to 1, .33 to .66 w2 = w draw mp left title "Pulmonary Artery Pressure (mmHg)" size .007 frame 0 to 1, 0 to .33 view 

The view command causes the graphs shown below to appear on the computer display:

/* read (time,pressure)-ordered pairs from press.dat */ pdat = read(press,71,2)      /* format = 71 rows, 2 columns */ /* define function for blood pressure as linear interpolation into pdat, and compute pressure for 1 cardiac cycle */ fct p(t) = lookup(pdat,t) ap = points(p,0:pdat[71,1]!100) /* define function for arterial flow under 2-element Windkessel model */ fct i(t) = p(t)/r+c*p't(t) /* set r,c constants for systemic arterial system and compute aortic flow for 1 cardiac cycles */ r = .9000                /* systemic peripheral resistance in (mmHg/cm^3/sec) */ c = 1.0666              /* systemic arterial compliance in (cm^3/mmHg) */ mi = points(i,0:pdat[71,1]!100) /* draw aortic pressure and aortic flow */ draw ap draw pdat lt none pt crosspt ptsize .005 left title "Aortic Pressure (mmHg)" size .01 no framebox frame 0 to 1, 0 to .5 w1 = w draw mi left title "Aortic Flow (ml/sec)" size .01 frame 0 to 1, .5 to 1. no framebox view 

The following graph results:

/* read digitized from MOL.DAT: time in seconds in column 1, pressure in mmHg in column 2, flow in ml/sec in column 3. */ dt = read(mol,65,3) n = nrows(dt) tv = dt col 1                         /* time in seconds */ pdat = tv '(dt col 2)     /* time,pressure ordered pairs */ idat = tv '(dt col 3)         /* time,flow ordered pairs */ /* define flow function and evaluate for one cardiac cycle */ fct i1(t) = lookup(idat,mod(t,tv[n])) af = points(i1,0:tv[n]!100) /* define function for evaluating Akaike criterion given s=sum-of-squares, nd=number of data points, and np=number of parameters in model */ fct aic(s,nd,np) = nd*log(s)+2*np /* 2-element Windkessel model with rw = peripheral resistance and cw = arterial compliance. Calculate values for parameters following Molino et al: peripheral resistance = mean pressure/mean flow; arterial compliance = time constant/peripheral resistance, where time constant = .529 seconds. */ fct pw't(t) = (i1(t)-pw(t)/rw)/cw init pw(0) = pdat[1,2] rw = mean(pdat col 2)/mean(idat col 2) cw = .529/rw mw2 = points(pw,0:tv[n]!100) maxiter = 0 fit (cw), pw to pdat a[1] = aic(sosq,n,2) maxiter = 20 /* 3-element Windkessel model with rb1 = valve resistance, rb2 = peripheral resistance, and cb = arterial compliance */ fct pb't(t) = rb1*i1't(t)+((1+rb1/rb2)*i1(t)-pb(t)/rb2)/cb init pb(0) = pdat[1,2] rb1 = 0 rb2 = rw cb = cw fit (rb1), pb to pdat a[2] = aic(sosq,n,3) mb3 = points(pb,0:tv[n]!100) /* 4-element Windkessel model with lm = inertial mass of fluid, rm1 = valve resistance, rm2 = peripheral resistance, and cm = arterial compliance. */ fct pm't(t)=lm*i1''t(t)+(rm1+lm/(cm*rm2))*i1't(t)+((1+rm1/rm2)*i1(t)-pm(t)/rm2)/cm init pm(0) = pdat[1,2] rm1 = rb1 rm2 = rb2 cm = cb lm = 0          /* in grams/cm^4 */ constraints q1 = {lm > 0} fit (lm), pm to pdat constraints q1 a[3] = aic(sosq,n,4) mm4 = points(pm,0:tv[n]!100) /* draw graphs */ draw af draw idat lt none pt circle ptsize .005 left title "Flow (ml/sec)" size .01 frame 0 to 1, .66 to .99 no framebox w1 = w draw mw2 draw mw2 row 1:100:15 lt none pt triangle ptsize .01 draw mb3 draw mb3 row 1:100:15 lt none pt square ptsize .01 draw mm4 draw mm4 row 1:100:15 lt none pt star ptsize .01 draw pdat draw pdat row 1:n:2 lt none pt crosspt ptsize .01 left title "Pressure (mmHg)" size .01 title "Pressure data = +" at (.2,.85) ffract size .01 title "2-element = '27TC'R, AIC = "+a[1] at (.2,.8) ffract size .01 title "3-element = '27TB'R, AIC = "+a[2] at (.2,.75) ffract size .01 title "4-element = '27TE'R, AIC = "+a[3] at (.2,.7) ffract size .01 frame 0 to 1, 0 to .66 no framebox view plot in molwtp.ps 

The resulting graph is: