Applications of Biomedical Engineering in Anatomy, Physiology, Medicine & Hospital Management

1. Applications of Biomedical Engineering in Anatomy, Physiology, Medicine & Hospital Management by Dhanjoo N. Ghista UNSW Asia, Singapore • Anatomy: Spinal Disc as an intrinsically optimal structure • Non-dimensional Physiological Index • Cardiac Fitness based on Treadmill Heart Rate Variation • Medical Physiology: A Non-dimensional diabetes index with respect to Oral- • Glucose-tolerance Testing • Cardiology: LV contractility based on its normalized wall stress • ICU: Indicator for lung-status in mechanically ventilated COPD patients • (using lung ventilation modeling and assessment) • Monitoring: Noninvasive determination of Aortic pressure, aortic stiffness and • peripheral resistance • Coronary Bypass Surgery: candidacy • Hospital Management

2. §Geometry and Deformation Variables of the Spinal Disc

3. Vertebral end-plate F -z σf b z a h Normal stresses f & z equilibrating the applied force P

4. Analysis for Vertical Loading For a vertically applied force P , (1) the axial stresses in the annulus. where is the hydrostatic pressure in the NP fluid, and We put (2)

5. We obtain (3) NP fluid pressure, (4) (5) in the annulus (6)

6. Axial deformation, Radial deformation The disc deformations are obtained as: (7) (8) wherein Ec = E - Eo = kσ. (9) Now, as the magnitude of the load P increases, the pressure p in nucleus-pulposus fluid also increases. Then, as p increases, so does the modulus Ec according to eqn (9)

7. = constant (10) Based on eqn (10), we have from eqn(s) (7 & 8), ,a constant and ,a constant This means that, irrespective of increase in the value of load P, the disc deformations and remain constant, and only depend on the constitutive property parameter k. This is the novelty of the intrinsic design of the spinal disc!

8. § Clinical Evaluation of Physiological Systems in terms of Non-dimensional Physiological Indices Non-dimensional numbers (made up of several phenomenon related terms) are employed to characterise – disturbance phenomena. For example, in a cardiovascular fluid-flow regime, the Reynold’s number (1) (V : flow velocity, D: diameter, : fluid viscosity, : fluid density) characterizes turbulent flow, which can occur in the ascending aorta when the aortic valve is stenotic (giving rise to murmurs) and accentuated in the case of anaemia (decreased blood viscosity).

9. Integration of a number of isolated but related events into one non-dimensional physiological index (NDPI) can help to characterise an abnormal state of a particular physiological system.

10. §Assessing Cardiac Fitness based on Treadmill HR Variation

11. in this procedure, the cardiac fitness model consists of a first-order differential-equation system model describing the heart rate (HR) response (y) to exertion (exercise, jogging etc) monitored in terms of the work-load, where y is defined as follows : (1) The subject is exercised on the treadmill for a period of time te (minutes) at a constant work load (W), while the HR(t) (and hence y) is monitored. Now we develop a model to simulate (i) the y(t) response during exercise, i.e. during t ≤ te, and (ii) thereafter for y(t) decay, after the termination of exercise. In a way, te represents the exercise endurance of the subject.

12. For a person, the model equation for HR response is represented by: (2) For t≤te, the solution is given by: (3) For the recovery period (t ≥ te ) , the solution of eqn. (2) is : (4) where k1 and k2, are the model parameters, which can serve as cardiac-fitness parameters (in min-1).

13. Non-dimensional Cardiac Fitness Index A typical y(t) response is illustrated in Fig 1. The parameters k1 and k2 can be continued into a single nondimensional cardiac-fitness index (CFI): (5) According to this formulation of CFI, a healthier subject has (i) greater k1 (i.e., slower rate-of-increase of HR during exercise (ii) greater k2 (i.e., faster rate-of-decrease of HR after exercise) (iii) greater te (i.e., exercise endurance), and hence (iv) higher value of CFI.

15. Now, we need to evaluate CFI for a big spectrum of patients, and then compute its distribution curve, to determine the efficiency of this index, in order to yield distinct separation of CFI ranges for healthy subjects and unfit patients. This CFI can also be employed to assess improvement in cardiac fitness following cardiac rehabilitation regime. This CFI is non-dimensional, and it can be useful to clinicians as they are able to predict the heart condition or fitness performance of a person by referring to the value of a single index value.

16. §MEDICAL PHYSIOLOGY : A NON-DIMENSIONAL DIABETES INDEX WITH RESPECT TO ORAL-GLUCOSE-TOLERANCE TESTING

17. For Oral-glucose Tolerance Test simulation (entailing digestive & blood-pool chambers), the differential equation is as follows : wherein ( = ) is the natural oscillation-frequency of the system, A is the attenuation or damping constant of the system,  = is the (angular) frequency of damped oscillation of the system, = 2A/Td = , with (λy) representing the proportional-control term of blood-glucose concentration (y). (λTdy’ )is the derivative feedback control term with derivative-time of Td represents the injected glucose bolus (1)

18. The input to this system is taken to be the impulse function due to the orally ingested glucose bolus [G], while the output of the model is the blood-glucose concentration response y(t). For an impulse glucose-input, a normal patient’s blood-glucose concentration data is depicted in Figure 1 by open circles. Based on the nature of this data, we can simulate it by means of the solution of the Oral-glucose regulatory (second-order system) model, as an under-damped glucose-concentration response curve, given by: y(t) = (G/ω)e-Atsin ωt, (2) wherein A is the attenuation constant, is the damped frequency of the system, the natural frequency of the system = ωn, and λ = 2A/Td. The model parameters λ and Td are obtained by matching eqn.(1) to the monitored glucose concentration y(t) data (represented by the open circles). The computed values of parameters are: λ = 2.6hr-2, Td = 1.08hr. This computed response is represented in Figure 1 by the bottom curve, fitting the open-circles clinical data.

19. Parametric Identification (sample calculation for Normal Test Subject No.5) y (1/2) = (G/)e – A/2 sin /2 = 0.34 g/L y (1) = (G/)e – A sin  = 0.24 g/L y (2) = (G/)e – 2Asin 2 = -0.09 g/L Using trignometry relations, we get A = 0.8287 hr –1  = n2 = A2 + 2 = (0.82875)2 + (2.0146)2 = 4.7455 hr– 2, Td = 2A/ = 0.3492 hr Upon substituting the above values of  and Td , the value of the third parameter, G = 1.2262 g (l)-1 hr-1

20. For a diabetic subject, the blood-glucose concentration data is depicted by closed circles in Fig 1. For the model to simulate this data, we adopt the solution of model eqn(17), as an over-damped response function: y(t) = (G/ω) e-At sinh ωt (3)

21. The soultion (y =(G/ ω)e-At sinhωt) is made to match the clinical data depicted by closed circles, and the values of λ and Td are computed to be 0.27hr-2 & 6.08hr, respectively. The top curve in Figure 1 represents the blood-glucose response curve for this potentially diabetic subject. The values of Td , λ and A for both normal and diabetic patients are indicated in the figure, to provide a measure of difference in the parameter values. hr-2 hr gm lit-1 hr-1 hr-2 hr gm lit-1 hr-1

22. It was found from these calculations that not all of the normal test subjects’ clinical data could be simulated as under-damped response. Similarly, not all the diabetic test subjects’ clinical data corresponded to over-damped response. However it was found that the clinical data of these test subjects (both normal and diabetic) could indeed be fitted by a critically-damped glucose-response solution of the governing equation. y(t) = G t e – At (4) for which,  = 0, n2 = A2 =  , and Td = 2A/ = 2

23. Clinically-based Diagnosis: The blood glucose ‘normal’ values, used for the clinical studies, were: Fasting : 70 to 115 mg/dl , At 30th min. : less than 200 mg/dl , At 1st hour : less than 200 mg/dl, At 2nd hour : less than 140 mg/dl, • Modeling-based Diagnosis: The test subjects have been classified into four categories: Normal-test subjects based on under-damped model-response; Normal test-subjects based on critically-damped model-response; Diabetic test-subjects based on critically-damped model-response; Diabetic test-subjects based on over-damped model response;

24. NON-DIMENSIONAL NUMBER FOR DIAGNOSIS OF DIABETES We decided to develop a unique diabetes index number (DIN) to facilitate differential diagnosis of normal and diabetic states as well as diagnose supposedly normal but high (diabetic) risk patients and diabetic patients in early stages of the disorder. DIN = (5) wherein, y(max) = maximum blood glucose value in gram/liter G = glucose dose administered to the system in gram/liter hour A = attenuation constant in 1/hour Td = derivative-time (+) in hours T(max) = the time at which y(max) is attained in hour This non-dimensional number DIN consists of the model parameters (A & Td) or (A & ωn) or (λ & Td)

25. The DIN values for all four categories were computed from equation (5). A distribution plot of the DIN is plotted, wherein the DIN is classified into sections with 0.2 increments ( for all the four categories of subjects ) and the number of subjects which fall into these sections (frequency) is determined. In the distribution plot (shown in Fig 2), the DIN values 0-0.2 is designated as range 1, the DIN 0.2-0.4 is range 2, 0.4-0.6 is range 3, and so on up to DIN 2.2-2.4, which is range 12. Fig. 2 DIN distribution plot of all the four categories subjects

26. As can be seen from the figure, normal (i.e., non-diabetic) subjects with no risk of becoming diabetic, will have DIN value less than 0.4, or be in the 1 – 2 range. Distinctly diabetic subjects will have DIN value greater than 1.2, or be in the 7 – 12 range categories. Supposedly, clinically-identified normal subjects who have DIN values between 0.6 and 1.0, or are in the 3 – 5 range, are at risk of becoming diabetic. On the other hand, clinically-identified diabetic subjects with DIN value between 0.6 – 1.2, or in the 3 – 6 range category are border-line diabetics, who can become normal (with diet control and treatment). Fig. 2 DIN distribution plot of all the four categories subjects

27. §Cardiology: LV Contractility Index based on Normalized Wall-Stress The traditional cardiac contractility index (dP/dt)max requires cardiac catheterization

28. Since LV Pressure is developed by LV wall stress (based on sarcomere contraction), we develop a contractility index based on  (normalized with respect to LV pressure) For a thick walled spherical model of LV, the differential stress: (1)

29. (2) (3)

30. Fig 1: Pressure and volume case study data

31. For the data shown in the figure 1 For a non-dimensional index,

32. Comparison of new & traditional indices

33. § ICU : INDICATOR FOR LUNG-STATUS IN MECHANICALLY VENTILATED COPD PATIENTS (USING LUNG VENTILATION MODELING AND ASSESSMENT) In chronic-obstructive-pulmonary-disease (COPD), elevated airway resistance and decreased lung compliance (i.e. stiffer lung) make breathing difficult. After these patients are mechanically ventilated, there is a need for accurate predictive indicators of lung-status improvement, for ventilator discontinuation through stepwise reduction in mechanical support, as and when patients are increasingly able to support their own breathing, followed by trials of unassisted breathing preceding extubation, and ending with extubation.

34. Governing Equation for V respect to PN(t) V.Pal Lung-Alveolar Model • wherein: • the driving pressure = PL; Pe= the end-expiratory pressure; • net pressure PN = PL– Pe • the parameters of the governing equation (1) are lung compliance (C) and airflow- • resistance (R), with both R & C being instantaneous values • V= V(t) – Ve (wherein Ve is the end expiratory lung volume)

35. V t t Pk PN Ptm B Pe t Figure 1: Lung ventilatory model data shows air-flow (V) and volume (V ) and net pressure (PN). Pause pressure (Ptm )occurs at tm, at which the volume is maximum (TV = tidal volume). Δt is the phase difference between the time of maximum volume and peak pressure (Pk); it also the time lag between the peak and pause pressures. B is the amplitude of the net pressure waveform PN applied by the ventilator. Thus, PN oscillates about Pe ,with amplitude B.

36. V t t Pk PN Ptm B Pe t We measure Peak pressure (Pk), Pause pressure (Ptm), tm & (or tm ) Then

37. Let B be the amplitude of the net pressure wave form applied by the ventilator, Cabe the average dynamic lung compliance, Ra the average dynamic resistance to airflow, the driving pressure PL =Pe+ B sin (ωt), and the net pressure PN = B sin (ωt). The governing equation (1) then becomes: The volume response to PN, the solution to eqn (2), is given by: • wherein: • (i)ka (=RaCa)is the average time constant, • the integration constant H is determined from the initial conditions, • the model parameters are Ca and ka(i.e. Ca and Ra), and • (iv)ω is the frequency of the oscillating pressure profile applied by the ventilator

38. An essential condition is that the flow-rate is zero at the beginning of inspiration and end of expiration. Hence, applying this initial condition of dV/dt =0 at t=0 to our differential eqn (3), the constant H is obtained as: Then from eqn (3) & (4), we obtain:

39. From equation (6), we obtain (by neglecting the following expression for ka: Evaluating parameters Ra & Ca : For evaluating the parameter ka (RaCa) , we will determine the time at which V (t) is maximum and equal to the tidal volume (TV), Hence putting dV/dt =0 in eqn (5), we obtain:

40. From eqn (5 & 6): At t=tm, the second term, Hence, eqn (8) becomes:

41. In eqn (10), if we put Then based on eqns (6 & 7), Then, based on equation (12), equation (10) becomes:

42. ( ) 16 Determining Lung-Compliance (Ca) and Air-Flow Resistance (Ra) From (13 & 7), Hence, from eqns (7 & 13), the average value of airflow-resistance (Ra) is: For our patients, the computed ranges of the parameters are: Now, that we have determined the expressions for Raand Ca, the next step is to develop an integrated index incorporating these parameters.

43. Formulating a Lung Ventilatory Index (LVI) incorporating Ra and Ca Let us obtain an order-of-magnitude values of LVI, for a mechanically ventilated COPD patient in acute respiratory failure: wherein,

44. Now, let us obtain an order-of-magnitude of LVI (by using representative computed values of Ra, Ca, RF, TV, and Pk) as abovefor a COPD patient with improving lung-status just before successful discontinuation. wherein, Hence, for LVI to reflect lung-status improvement in a mechanically ventilated COPD patient in acute respiratory failure, it has to decrease to the range of LVI for an outpatient COPD patient at the time of discontinuation.

45. Successful discontinuation cases Failed discontinuation cases Figure 2: Distribution of LVI at discontinuation for patients with failed and successful discontinuation. For the 6 successfully-discontinued cases, the LVI was (2900) ± (567) (cmH2O/L)3; for the 7 failed-discontinuation cases the LVI was (11400) ± (1433) (cmH2O/L)3. It is observed that LVI enables clear separation between failed and successful discontinuation.

46. § MONITORING: Noninvasive determination of Aortic pressure, Aortic modulus (stiffness) and peripheral resistance) The aortic blood pressure waveform contains a lot of information on how the LV contraction couples with the aortic compliance and peripheral resistance. Since accurate measurement of aortic blood pressure waveform requires catheterization of the aorta, we have developed a noninvasive method to determine the aortic pressure profile along with the aortic volume-elasticity and peripheral resistance.

47. The input to the model consists of auscultatory cuff diastolic and systolic pressures, along with the cineangiographically measured ejection volume-time profile (or volume input into the aorta). The governing differential equation for pressure response to LV outflow rate I(t) into the aorta is given by (1) where (i) ma is the aortic volume elasticity (dP/dV), (ii) Rp is the resistance to aortic flow (=P/Q), (iii) =ma/Rp , (iv) I(t)is the monitored inflow rate into the aorta, and (v) Ta is the flow-acceleration period. This governing equation is solved for measured I(t) and dI/dt during the systolic phase from time T1 to T3 (Fig. 1). For the diastolic-phase solution from time T3 to T4, the right-hand side is zero.

48. Solution equations: • During diastolic phase (2) At t=T, Pd(T)=auscultatory Pad =P1 • During systolic phase (3) where Also, at t=T2, dPs/dt=0, at t=T2, Ps(T2)=auscultatory Pas=P2; at t=T3, Ps(t=T3)=Pd(t=T3)

49. Equations to be solved: ● (4) ● At t=T2, Ps=P2. Hence, from equations (3 & 4), we get: (5) At t=T3, Ps(t=T3)=Pd(t=T3) ● (6) Now solve equations (4, 5 & 6), to determine the unknown parameters ma, Rp (and T2)

50. We have determined the expressions for the aortic pressure during the systolic and diastolic phases, by solving the governing equation (1), for the monitored LV outflow rate (or input into the aorta) I(t), using (i) the monitored auscultatory diastolic pressures (Pad), to serve as the boundary condition at the beginning of the systolic-phase solution (at time T1) and at the end of the diastolic-phase solution (at time T4), (ii) the monitored auscultatory systolic pressure (Pas =P2) to represent the maximum value of the systolic-phase solution. Because the pressure solution of equation (1) is a function of ma and Rp, we first determine the values of these parameters by making the solution satisfy the above-mentioned boundary conditions, which in turn yields the pressure profile.