BLOOD CIRCULATION AND AQUEOUS HUMOR FLOW IN THE EYE: MULTI-SCALE MODELING AND CLINICAL APPLICATIONS A Dissertation Submitted to the Faculty of Purdue University by Simone Cassani In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy August 2016 Purdue University Indianapolis, Indiana ii To Lucia, Dad, Mom, Alessandra and Martina. iii ACKNOWLEDGMENTS First, I would like to express my gratitude to my advisor Professor Giovanna Guidoboni for her help and guidance during this PhD on both the academic and human levels. Working with a mentor of your knowledge and experience has been invaluable to me. You taught me that, with determination and hard work, every objective can be reached and I owe you the most for the researcher I am today. I also would like to devote a special thank you to my co-advisors Professors Julia Arciero and Alon Harris. Your expertise and your help played an important role in my PhD research and I am grateful I had the possibility to work with you. My sincere thanks also go to Professors Christophe Prud’homme, Marcela Szopos and Riccardo Sacco and Doctors Aurelio Mauri and Brent Siesky. Our collaborations in Strasbourg, Milan and Indianapolis allowed me to widen my research and enrich my skills-set. I want to thank my family. You have always supported me and I am proud of you. I want to thank you for the sacrifices you made for me and I want you to know that I couldn’t have reached this achievement without you. I thank all my friends (on both sides of the Atlantic Ocean) for all the moments we spent together at school, on skype, at home and in bars. I know I am very lucky to be surrounded by wonderful people like you. Last, I would like to thank my wife-to-be Lucia. You have always been there to celebrate my accomplishments and to encourage me when all I could see was a gray sky. We started this journey together and after many sleepless nights, some struggles and a lot of work we finished it together. I hope I made you proud as much as I am proud of you. Having you by my side is the greatest gift I could have ever asked for. (h). iv TABLE OF CONTENTS Page LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 The physiology of the eye . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Mathematical modeling of the retinal circulation . . . . . . . . . . . . . . 7 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . 9 2.1.2 Control state . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.3 Time-profiles of the input and output pressure waves . . . . 16 2.1.4 IOP-induced stress in the lamina cribrosa . . . . . . . . . . . 17 2.1.5 Variable passive resistance . . . . . . . . . . . . . . . . . . . 19 2.1.6 Active variable resistance . . . . . . . . . . . . . . . . . . . . 29 2.2 Solution algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3.1 Model validation . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3.2 Theoretical investigations . . . . . . . . . . . . . . . . . . . 35 2.3.3 Comparison between model predictions and clinical data . . 38 2.3.4 Trabeculectomy simulations . . . . . . . . . . . . . . . . . . 41 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.4.1 Model validation . . . . . . . . . . . . . . . . . . . . . . . . 46 2.4.2 Theoretical investigations . . . . . . . . . . . . . . . . . . . 47 2.1 2.4 v Page 2.4.3 Comparison between model predictions and clinical data . . 49 2.4.4 Trabeculectomy simulations . . . . . . . . . . . . . . . . . . 50 3 Coupled model for vascular regulation in the retina . . . . . . . . . . . . 52 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.1.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . 53 3.1.2 Modeling of vascular regulation . . . . . . . . . . . . . . . . 56 3.1.3 Control state . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.2 Solution algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.3.1 Effect of model coupling . . . . . . . . . . . . . . . . . . . . 65 3.3.2 Comparison with experimental measurements for changes in MAP and IOP . . . . . . . . . . . . . . . . . . . . . . . . . 67 Comparison with experimental measurements for changes in oxygen demand . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.1 3.3.3 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.4.1 Effect of model coupling . . . . . . . . . . . . . . . . . . . . 74 3.4.2 Comparison with experimental measurements for changes in MAP and IOP . . . . . . . . . . . . . . . . . . . . . . . . . 75 Comparison with experimental measurements for changes in oxygen demand . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.4.3 4 Modeling oxygenation within multiple retinal layers . . . . . . . . . . . . 78 Oxygen transport in the retinal tissue . . . . . . . . . . . . . . . . . 78 Model reduction . . . . . . . . . . . . . . . . . . . . . . . . . 81 Oxygen transport in the retinal microvasculature . . . . . . . . . . 83 Model reduction . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.3 Coupling between retinal tissue and retinal microvasculature . . . . 95 4.4 Solution algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.5 Preliminary results and discussion . . . . . . . . . . . . . . . . . . . 102 5 Mathematical modeling of aqueous humor flow . . . . . . . . . . . . . . . 105 4.1 4.1.1 4.2 4.2.1 vi Page 5.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 VITA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 vii LIST OF TABLES Table Page 2.1 Geometric and physical parameters for the CRA and CRV. . . . . . . . 12 2.2 Computed values of the vascular resistances at the control state. . . . . 13 2.3 Geometric and physical parameters used in the computation of control state resistances for the arterioles, capillaries and venules compartments in the control state. Adapted from [59]. . . . . . . . . . . . . . . . . . . 14 2.4 Computed values of intraluminal blood pressure at the control state. . 15 2.5 Computed values of vascular volume and capacitance at the control state. 16 2.6 Radius and thickness values used in the mathematical model for the lamina cribrosa and sclera . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Model values of Young’s modulus and shear modulus for the lamina cribrosa as a function of the effective stress . . . . . . . . . . . . . . . . . . . . 18 2.8 Representative cases for model simulations . . . . . . . . . . . . . . . . 34 2.9 Changes in peak systolic velocity, end diastolic velocity, resistivity index and mean flow velocity in the central retinal artery as IOP is reduced by trabeculectomy, deep sclerectomy and model simulations . . . . . . . . 40 2.7 2.10 Representative cases for trabeculectomy model simulations . . . . . . . 42 2.11 Changes in IOP and CRA hemodynamics following trabeculectomy [88]. 43 3.1 Computed values of the vascular resistances at the control state. . . . . 57 3.2 Values of the constant parameters for vascular regulation mechanisms. 59 3.3 Control state values for retinal vascular network compartments. . . . . 62 3.4 Computed values of intraluminal blood pressures at the control state. . 63 3.5 Summary of vascular response to flicker stimulation and light-dark adaptation studies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.1 Parameters for the oxygen transport in the retinal tissue . . . . . . . . 82 4.2 Parameters for the oxygen transport in the retinal microvasculature . . 86 4.3 Geometrical and physiological parameter of the network in Figure 4.2 . 95 viii Table 5.1 5.2 Page Control state values for the parameters in the model for aqueous humor flow in equation (5.10). . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Mean values, standard deviations and skewness of the distribution of IOP resulting from the sensitivity analysis of the mathematical model for four cases of clinical interest . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 ix LIST OF FIGURES Figure Page 1.1 Anatomy of the eye . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Optic nerve canal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Retinal blood flow response to changes in pressure and oxygen demand in the case of normal retinal vascular regulation . . . . . . . . . . . . . . 3 2.1 Network model for the retinal vasculature . . . . . . . . . . . . . . . . 8 2.2 Relationship between the DN model by Takahashi et al. [59] and the chosen geometry of our model . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Typical color Doppler imaging measurement of the blood velocity in the CRA and CRV; time profiles of the blood velocity in the CRA and CRV at the control state; and time profile of the inlet and outlet pressures at the control state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Schematic representation of geometry and boundary conditions of the elasticity problem for the lamina cribrosa . . . . . . . . . . . . . . . . . . . 19 2.5 Representative cylinder Ω in the Cartesian coordinate system (x, y, z) . 20 2.6 Radial displacement on cross-section Σ . . . . . . . . . . . . . . . . . . 25 2.7 Graph of tube law from experimental results for a collapsible tube. . . 26 2.8 Arteriolar resistance and normalized flow for model simulations with functional or absent autoregulation for various ocular perfusion pressures (OPP) 30 Schematic representation of the translaminar portion of the central retinal vessels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.10 Comparison of model predicted values of blood velocity, volumetric blood flow, and normalized total retinal blood flow with measured data . . . 35 2.11 Model predicted values of total retinal blood flow; peak systolic velocity, end diastolic velocity, and resistivity index in the central retinal artery as IOP varies between 15 and 45 mmHg for six different theoretical patients 36 2.12 Model predicted ranges of IOP in 15 − 30 mmHg for which the values of total retinal blood flow, PSV, EDV and RI are within ±3% of their control state values for six different theoretical patients . . . . . . . . . . . . . 37 2.3 2.4 2.9 x Figure Page 2.13 Model predicted values of blood pressure in CRA, arterioles, capillaries, venules, and CRV for IOP = 15, 30, and 45 mmHg for six different theoretical patients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.14 Comparison between PSV and EDV measured by Harris et al. [11] in the CRA and the model predicted values in the NBPwAR and NBPwoAR cases; and comparison between the percent change in CRA mean flow velocity measured by Findl et al. [9] and the model predicted values in the NBPwAR and NBPwoAR cases . . . . . . . . . . . . . . . . . . . . . . 41 2.15 Scatter plot of the velocity difference in the CRA for the patients in [88] 44 2.16 Model predictions and clinical data value of CRA velocity before and after trabeculectomy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.17 Model predictions of average retinal blood flow in a cardiac cycle for the four theoretical patients. . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.1 Schematic for retinal vascular network in the coupled model . . . . . . 54 3.2 Krogh cylinder model schematic representation . . . . . . . . . . . . . 60 3.3 Model predicted values of P1,2 and P4,5 , pressure drop along the CRA, and total resistance in the CRA for the microcirculation model and coupled model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Predicted values of retinal blood flow (normalized) and predicted values for venous oxygen saturation as MAP is varied (65 − 165 mmHg) using the microcirculation model and coupled model . . . . . . . . . . . . . . 67 Model predicted values of retinal blood flow (normalized) as MAP is varied (20 − 165 mmHg) using the coupled model. Vascular regulation mechanisms are turned on and off to assess their individual contributions . . 68 Normalized values of blood flow predicted by the coupled model and microcirculation model are compared with data from clinical studies [39,120,121] 69 Coupled model and microcirculation model predictions of normalized blood flow as IOP is varied between 7.5 − 31 mmHg are compared with experimental measures [39, 120] at low-, normal- and elevated- blood pressure values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Large arteriole blood flow vs. oxygen demand. The model predictions for different vascular regulation mechanisms are compared with clinical data [30] at baseline and during flicker stimulation . . . . . . . . . . . . 72 4.1 Tissue layers of the retina . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.2 Retinal microvasculature . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.4 3.5 3.6 3.7 3.8 xi Figure Page 4.3 Representative cylinder Ω in cylindrical coordinates (ζ, θ, s) . . . . . . 85 4.4 Schematic of the solution algorithm for the coupled problem of oxygen transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 Model predictions of blood oxygen partial pressure as a function of position in the network in the six microcirculation compartments at control state 103 Model predictions of blood oxygen saturation as a function of position in the network in the six microcirculation compartments at control state . 104 Model predictions of tissue oxygen partial pressure distribution in the six retinal layers in Figure 4.1. . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.1 Network model for production and drainage of aqueous humor. . . . . . 106 5.2 Probability density function of intraocular pressure and Sobol indices resulting from the sensitivity analysis performed on the mathematical model 115 4.5 4.6 4.7 xii ABBREVIATIONS ACA Anterior ciliary arteries AH Aqueous humor AR Autoregulation ATP Adenosine triphosphate BP Blood pressure cBP Blood pressure in the capillaries of the ciliary body C Capillaries CDI Color Doppler imaging CO2 Carbon dioxide CRA Central retinal artery CRV Central retinal vein DN Dichotomous network DP Diastolic arterial blood pressure EDV End diastolic velocity EVP Episcleral venous pressure IOP Intraocular pressure LA Large arterioles LV Large venules LC Lamina cribrosa MAP Mean arterial pressure at the level of the brachial artery MR Metabolic regulation MFV Mean flow velocity O2 Oxygen OA Ophthalmic artery xiii OPP Ocular perfusion pressure PCA Posterior ciliary arteries PSV Peak systolic velocity RBC Red blood cell RI Resistive index RLTp Retrolaminar tissue pressure SA Small arterioles SP Systolic arterial blood pressure SV Small venules TB Trabeculectomy VR Vascular regulation xiv ABSTRACT Cassani, Simone PhD, Purdue University, August 2016. Blood Circulation and Aqueous Humor Flow in the Eye: Multi-Scale Modeling and Clinical Applications. Major Professor: Giovanna Guidoboni. Glaucoma is a multi-factorial ocular disease associated with death of retinal ganglion cells and irreversible vision loss. Many risk factors contribute to glaucomatous damage, including elevated intraocular pressure (IOP), age, genetics, and other diseases such as diabetes and systemic hypertension. Interestingly, alterations in retinal hemodynamics have also been associated with glaucoma. A better understanding of the factors that contribute to these hemodynamic alterations could lead to improved and more appropriate clinical approaches to manage and hopefully treat glaucoma patients. In this thesis, we develop several mathematical models aimed at describing ocular hemodynamics and oxygenation in health and disease. Precisely we describe: (i) a time-dependent mathematical model for the retinal circulation that includes macrocirculation, microcirculation, phenomenological vascular regulation, and the mechanical effect of IOP on the retinal vasculature; (ii) a steady-state mathematical model for the retinal circulation that includes macrocirculation, microcirculation, mechanistic vascular regulation, the effect of IOP on the central retinal artery and central retinal vein, and the transport of oxygen in the retinal tissue using a Krogh cylinder type model; (iii) a steady-state mathematical model for the transport of oxygen in the retinal microcirculation and tissue based on a realistic retinal anatomy; and (iv) a steady-state mathematical model for the production and drainage of aqueous humor (AH). The main objective of this work is to study the relationship between IOP, systemic blood pressure, and the functionality of vascular autoregulation; the transport xv and exchange of oxygen in the retinal vasculature and tissue; and the production and drainage of AH, that contributes to the level of IOP. The models developed in this thesis predict that (i) the autoregulation plateau occurs for different values of IOP in hypertensive and normotensive patients. Thus, the level of blood pressure and functionality of autoregulation affect the changes in retinal hemodynamics caused by IOP and might explain the inconsistent outcomes of clinical studies; (ii) the metabolic and carbon dioxide mechanisms play a major role in the vascular regulation of the retina. Thus, the impairment of either of these mechanisms could cause ischemic damage to the retinal tissue; (iii) the multi-layer description of transport of oxygen in the retinal tissue accounts for the effect of the inner and outer retina, thereby improving the predictive ability of the model; (iv) a greater reduction in IOP is obtained if topical medications target AH production rather that AH drainage and if IOP-lowering medications are administrated to patients that exhibit a high initial level of IOP. Thus, the effectiveness of IOP-lowering medications depend on a patient’s value of IOP. In conclusion, the results of this thesis demonstrate that the insight provided by mathematical modeling alongside clinical studies can improve the understanding of diseases and potentially contribute to the clinical development of new treatments. 1 1. INTRODUCTION 1.1 The physiology of the eye The eye is a complex structure composed of different parts, as depicted in Fig- ure 1.1, that work together to achieve vision. Anterior chamber (aqueous humour) Posterior chamber Suspensory ligament of lens Cornea Pupil Uvea Iris Ciliary body Lens Choroid Sclera Vitreous humour Retinal blood vessels Hyaloid canal Retina Macula Fovea Optic nerve Optic disc Figure 1.1. Anatomy of the eye. By Rhcastilhos Schematic_diagram_of_the_human_eye_with_English_annotations.svg, Public Domain, via Wikimedia Commons. The sclera is the outermost part of the eye shell that protects the inner layers of the eye. The choroid is the mid layer of the eye shell, located between the sclera and the retina, that nourishes the outer part of the retina. The retina is the inner layer of the eye shell that absorbs the light and sends the visual stimulus to the brain; it 2 contains the macula and the fovea that are the regions responsible for visual acuity. The ciliary body is located in the anterior part of the eye; it contains the ciliary muscle and produces the aqueous humor (AH), a watery fluid that fills the space in front of the lens, maintains the pressure inside the eye (intraocular pressure, IOP), keeps the eye inflated and transports nutrients to the avascular regions of the eye. The vitreous is a gel-like substance that fills the region behind the lens. The main source of blood to the eye is the ophthalmic artery (OA), a branch of the inner carotid artery. The major branches of the OA are the anterior ciliary arteries (ACA), that nourish the ciliary body and the choroid, the central retinal artery (CRA), that nourishes the inner part of retina, and the posterior ciliary arteries (PCA), that nourish the choroid. Lamina cribrosa Retina Choroid Sclera Posterior short ciliary artery and vein Bundles of optic nerve Central retinal artery and central retinal vein Figure 1.2. Optic nerve canal. The central retinal artery (CRA) and the central retinal vein (CRV) run into the optic nerve canal piercing through a collagen structure called lamina cribrosa (LC). Adapted from: by Henry Vandyke Carter, Public Domain, via Wikimedia Commons, from [1], Plate 880. 3 The CRA runs through the optic nerve canal, as depicted in Figure 1.2, piercing a collagen structure called the lamina cribrosa (LC) that is responsible for maintaining the pressure difference between the IOP (inside of the eye) and the retrolaminar tissue pressure (in the optic nerve canal). When the CRA reaches the inside of the eyeball it branches into 4 main retinal arteries. These vessels branch into smaller arteries, arterioles, and capillaries that supply blood and nutrients to the retinal tissue. Blood exits the capillaries and drains into venules, small veins, and the 4 major retinal veins. These vessels converge into the central retinal vein (CRV) that runs through the optic nerve canal parallel to the CRA, piercing the LC and converges into the ophthalmic retinal blood blood flow retinal retinal flow retinal blood blood flow vein and cavernous sinus. pressure pressure oxygen oxygenconsumption demand Figure 1.3. Retinal blood flow response to changes in pressure (left) and oxygen demand (right) in the case of normal retinal vascular regulation. Under healthy conditions, the retina exhibits vascular regulation, which is the ability to adjust perfusion in response to alterations in blood pressure or tissue oxygen demand, as schematically represented in Figure 1.3. This change in perfusion is achieved through vascular responses (eg. myogenic, shear stress, metabolic and carbon dioxide responses) to stimuli that affect the vascular tone of the smooth muscle cells inside the wall of retinal arteries causing constriction or dilation of these walls. 4 1.2 Motivation Alterations in retinal and retrobulbar hemodynamics are associated with many ocular and systemic diseases, including glaucoma [2–4], age-related macular degeneration [5, 6] and diabetes [6, 7]. Furthermore, an elevated level of IOP has been established as a risk factor for the prevalence, incidence and development of glaucoma, a disease that is the second leading cause of blindness worldwide. Since the retinal vasculature is in direct contact with IOP, it seems very likely that IOP impacts retinal blood flow and velocity. However, clinical evidence of such an impact is not unanimous. Several clinical studies have shown a decrease in retinal and retrobulbar blood flow and velocity when IOP was increased [8–15], while other studies did not find any significant changes in retinal and retrobulbar hemodynamics with changes in IOP [16–28]. In the clinic, elevated IOP is treated with topical medications and/or with surgical procedures, such as trabeculectomy. IOP-lowering topical medications affect the production and/or drainage of AH, and may alter retrobulbar hemodynamics. For example, studies found slightly increased retrobulbar circulation in patients on these medications [8,12,13], while most of the studies found that treatment with topical medications did not have a significant effect on retrobulbar hemodynamics [16, 17, 20–24, 26–28]. During trabeculectomy surgery, IOP is reduced by removing part of the eye trabecular meshwork and adjacent structures, thereby increasing AH outflow. Most of the hemodynamic measurements obtained after trabeculectomy provide evidence of increased blood velocity as a result of lowering IOP [10,14,15]. However, other studies did not report any significant change in blood velocity in the central retinal artery after trabeculectomy [18, 25]. Only a few studies have been conducted to evaluate the effect of IOP elevation on ocular hemodynamics, and these have also provided variable results [9, 11, 19]. 5 It is hypothesized that in a disease state, some of the vascular regulation mechanisms might be impaired, thereby compromising the oxygenation in the retina. However, there is inconsistency in the scientific literature regarding the vascular response to changes in oxygen demand, as observed in flicker stimulation studies and light-dark adaptation studies [29–36]. These inconsistent clinical observations are likely due to the numerous factors, including arterial blood pressure [3, 37–41] and vascular regulation [3, 42–44], that influence the relationship between IOP and ocular hemodynamics, and the intrinsic difficulty of evaluating the individual contributions of these factors in a clinical setting. Mathematical modeling can be used to investigate the complex relationship among these factors and to interpret the outcomes of clinical studies. In this perspective, the aims of this thesis are: • To develop a reduced mathematical model for the retinal circulation that simultaneously accounts for (i) blood flow in the central retinal vessels; (ii) blood flow in the retinal microcirculation; (iii) phenomenological retinal blood flow regulation; (iv) biomechanical action of IOP on the retinal vasculature; and (v) time-dependent arterial blood pressure. The model will be used to investigate the relationship between IOP, blood pressure (BP) and vascular autoregulation (AR) and the impact of these factors on retinal hemodynamics (Chapter 2). • To characterize the difference between compressible and collapsible tubes that are used to describe the wall deformation in arteries and veins (Section 2.1.5). • To improve the description of vascular regulation in the model developed in Chapter 2 by coupling it with an existing mechanistic model for retinal vascular regulation (described in [45]). The model in [45] utilizes a Krogh-cylinder type model to account for level of oxygen demand in the tissue. The resulting model will include both the effect of IOP on the CRA and CRV, and the mechanistic description of vascular regulation (Chapter 3). 6 • To derive and implement a model for the transport of oxygen in the retinal microvasculature and retinal tissue based on a physiological representation of: (i) the retinal layers, (ii) the retinal vasculature, and (iii) the oxygen exchange between the vasculature and the retinal tissue (Chapter 4). • To investigate the processes that contribute to the level of IOP, and the uncertainties in the outcomes of administration of IOP-lowering topical medication using a model that describes the production and drainage of AH (Chapter 5). 7 2. MATHEMATICAL MODELING OF THE RETINAL CIRCULATION This Chapter introduces a mathematical model for the retinal circulation that simultaneously accounts for (i) blood flow in the central retinal vessels; (ii) blood flow in the retinal microcirculation; (iii) retinal blood flow autoregulation; (iv) the biomechanical action of IOP on the retinal vasculature; and (v) time-dependent arterial blood pressure. The model employs a hydraulic analogy to Ohm’s law in which blood flow in the central retinal vessels and retinal microvasculature is analogous to the current flowing through a network of resistances and capacitances. Variable resistances describe active and passive diameter changes due to blood flow autoregulation and IOP. Systemic arterial BP is an input to the model. Model outputs include total retinal blood flow and blood velocity in the CRA, CRV and microvasculature, which are compared with clinically-measurable quantities. In Section 2.1 the model and its governing equations are described; in Section 2.2 the algorithm to perform numerical simulations is defined; in Section 2.3 the mathematical model is used to understand the relationships among IOP, BP, autoregulation and retinal and retrobulbar hemodynamics, and to investigate the effect of trabeculectomy on retinal hemodynamics. In Section 2.4 the results are discussed and compared with clinical data. The material and results presented in this Chapter are partially published in [46, 47] 2.1 Methods The retinal circulation is described using the analogy between the flow of a fluid in a hydraulic network and the flow of current in an electric circuit, as depicted in Figure 2.1. The vasculature supplying the retina is divided into five main compart- 8 P1,2 venules arterioles capillaries P P3,4 R2a P2 R2b 2,3 R3a R3b R4a P4 R4b P4,5 C2 IOP P3 C4 R1d P5,post P1,post R1c P1,pre R1b C1 CRA intraocular region IOP LAMINA CRIBROSA P5,pre CRV C5 P1 R5a R5b R5c P5 R1a RLTp P1,in R5d RLTp translaminar region retrobulbar region P5,out Rin Rout Pin (t) Pout (t) Figure 2.1. Network model for the retinal vasculature. The vasculature is divided into five main compartments: the central retinal artery (CRA), arterioles, capillaries, venules, and the central retinal vein (CRV). Each compartment includes resistances (R) and capacitances (C). The intraocular segments are exposed to the IOP, the retrobulbar segments are exposed to the RLTp, and the translaminar segments are exposed to an external pressure that depends on the internal state of stress within the lamina cribrosa (gray shaded area). Diameters of venules and intraocular and translaminar segments of the CRA and CRV are assumed to vary passively with IOP, whereas arterioles are assumed to be vasoactive. ments: the CRA, arterioles, capillaries, venules, and the CRV. Using the analogy between hydraulic and electrical circuits, blood flow is modeled as current flowing through a network of resistors (R), representing the resistance to flow offered by blood vessels, and capacitors (C), representing the ability of blood vessels to deform and store blood volume. Resistors and capacitors have been labeled with numbers from 1 to 5 to distinguish between compartments, and additional alphabetic labels are used to distinguish between segments within the same compartment. For example, the CRA compartment is referenced by the label 1 and it includes a retrobulbar 9 segment with resistances R1a and R1b , a translaminar segment with resistance R1c , and an intraocular segment with resistance R1d . We remark that the resistance of the retrobulbar segment is split into R1a and R1b to allow the capacitance C1 to act on the mean pressure of the retrobulbar portion of the CRA (P1 ). The blood flow through the retinal vascular network is driven by Pin and Pout , which represent the BPs upstream of the CRA and downstream of the CRV, respectively. The vascular segments are exposed to various external pressures depending on their position in the network. The intraocular segments are exposed to the IOP, the retrobulbar segments are exposed to the retrolaminar tissue pressure (RLTp), and the translaminar segments are exposed to an external pressure that depends on the internal state of stress within the lamina cribrosa. The IOP-induced stress within the lamina cribrosa is computed using a nonlinear elastic model described in [48]. The resistances of the venules and intraocular and translaminar segments of the CRA and CRV are assumed to vary passively with IOP, as detailed in Section 2.1.5. Arteriolar resistance varies actively to maintain a relatively constant blood flow despite changes in the ocular perfusion pressure (OPP). Ocular perfusion pressure is defined as OPP= (2/3) MAP − IOP, where MAP is the mean arterial pressure at the level of the brachial artery, MAP= (2/3)DP +(1/3) SP, and DP and SP are diastolic and systolic arterial BPs, respectively. The arrows in Figure 2.1 indicate the resistors that vary passively or actively. 2.1.1 Governing equations Ohm’s law states that the flow Q through a resistor is directly proportional to the pressure drop ∆P across the resistor, with a proportionality constant that is equal to the reciprocal of the resistance, namely Q= ∆P . R (2.1) 10 The flow Q through a capacitor is directly proportional to the time derivative of the product between the pressure drop across the capacitor and the capacitance, namely Q= d (C∆P ) . dt (2.2) Here, capacitance is assumed constant so (2.2) becomes Q=C d∆P . dt (2.3) By Kirchoff’s law, the following relationship must hold at every node in the network: volume change = flow in − flow out. (2.4) The application of Kirchoff’s law to the retinal vascular network shown in Figure 2.1 leads to the following system of ordinary differential equations for the nodal pressures P1 , P2 , P4 , P5 :     d(P1 − RLTp)    C1  dt       d(P2 − IOP)    C2 dt   d(P4 − IOP)    C4   dt        C d(P5 − RLTp)   5 dt = P1 − P2 Pin − P1 − Rin + R1a R1b + R1c + R1d + R2a = P1 − P2 P2 − P4 − R1b + R1c + R1d + R2a R2b + R3a + R3b + R4a = P2 − P4 P4 − P 5 − R2b + R3a + R3b + R4a R4b + R5a + R5b + R5c = P 4 − P5 P5 − Pout − . R4b + R5a + R5b + R5c R5d + Rout (2.5) The inlet and outlet pressures Pin and Pout vary with time along a cardiac cycle and, consequently, the pressures calculated via equation (2.5) are time dependent. The remaining nodal pressures in the circuit of Figure 2.1 can be computed via Ohm’s law in the corresponding segments. 2.1.2 Control state A control state for the system is defined to represent typical conditions of a healthy eye. Control values of any given quantity will be indicated with an overline bar. 11 External and systemic pressures We assume that for a healthy individual IOP = 15 mmHg [49–51], RLTp = 7 mmHg [52, 53], SP = 120 mmHg, and DP = 80 mmHg. Retinal blood flow The control value Q of the retinal blood flow can be estimated by applying Poiseuille’s law to the CRA, since the CRA is the only vessel supplying blood to the retinal vascular network. In the hypotheses of laminar flow and cylindrical geom2 etry, Poiseuille’s law yields Q = πV cra,max Dcra /8, where Dcra is the CRA diameter, and V cra,max is the average value of the CRA centerline velocity over a cardiac cycle. Measurements from human patients are used to define the control state values for Dcra = 175 µm [54], and V cra,max = 5.67 cm/s [11, 54], which yield a physiological value of Q = 6.8178 · 10−4 mL/s [54–58]. Resistance According to Poiseuille’s law, the resistance of a vessel is R= 128µL , πD4 (2.6) where D is the vessel diameter, L is the vessel length and µ is blood viscosity. Using this law and the geometric and physical parameters for the CRA and CRV reported in Table 2.1, the control state values of the CRA and CRV resistances are computed and summarized in Table 2.2. Determining the control state values of arteriolar, capillary, and venular resistance is more complex, because these compartments include a hierarchy of numerous vessels of various diameter. We adapt the dichotomous network (DN) model for the retinal microcirculation proposed in [59] to describe the hierarchical architecture of arterioles, capillaries, and venules. The DN model includes 14 levels of arterioles, 1 level of capillaries, and 14 levels of venules; each level includes a specific number of identical 12 Table 2.1. Geometric and physical parameters for the CRA and CRV. CRA Source CRV Parameters Value Value Source Diameter D, µm 175 [54] 238 [59]* Length L, mm, total 10 [60] 10 [60] Segment a 4.4 [60] 1 [61] Segment b 4.4 [60] 0.2 [62, 63] Segment c 0.2 [62, 63] 4.4 [60] Segment d 1 [61] 4.4 [60] Blood viscosity µ, cP 3.0 [64, 65] 3.24 [59]** Wall Young’s modulus E, MPa 0.3 [65–67] 0.6 [68] Wall Poisson’s ratio ν 0.49 [66, 67] 0.49 *** Wall thickness h, µm 39.7 [65, 69] 10.7 [70] *CRA/CRV diameter ratio is assumed to be the same as that for the first generation arteries/veins in [59]. **Blood viscosity in the CRV is assumed to be 8% larger than that in the CRA (as in [59]). ***Poisson’s ratio of the CRV wall is assumed to be the same as that for the CRA wall. parallel vessels. In our model, we characterize these 29 vascular levels into arterial, capillary and venous compartments according to vessel size. Vessels with diameter less than 6.5 µm are defined as capillaries, as shown in Figure 2.2. Vessel number, diameter, length, and blood viscosity for the 29 levels in the DN model are reported in Table 2.3 and the corresponding values of the lumped resistances calculated in the control state are summarized in Table 2.2. It is important to note that the viscosity values used in the model are effective viscosity values that are based on an empirical 13 Table 2.2. Computed values of the vascular resistances at the control state. Value, Resistance Value, mmHg · s/mL Resistance mmHg · s/mL Rin 2.25 · 104 R3b 5.68 · 103 R1a 4.30 · 103 R4a 3.11 · 103 R1b 4.30 · 103 R4b 3.11 · 103 R1c 1.96 · 102 R5a 3.08 · 102 R1d 9.78 · 102 R5b 6.15 · 101 R2a 6.00 · 103 R5c 1.35 · 103 R2b 6.00 · 103 R5d 1.35 · 103 R3a 5.68 · 103 Rout 5.74 · 103 relationship and that depend on vessel diameter. In this way, the model takes into account the corpuscular nature of blood. arterioles Our model DN model 1art ... capillaries 12art13art14art cap. 14ven13ven venules ... 1ven diameter < 6.5µm Figure 2.2. Relationship between the DN model by Takahashi et al. [59] and the chosen geometry of our model. The 29 vascular levels of the DN model have been divided into three model compartments for arterioles, capillaries, and venules according to vessel size. Rin and Rout represent the ocular vessels upstream of the CRA and downstream of the CRV. Their control values are determined by the control value Q of the total 14 Table 2.3. Geometric and physical parameters used in the computation of control state resistances for the arterioles, capillaries and venules compartments in the control state. Adapted from [59]. Arterioles in [59] Level Venules in [59] Number Diameter Length Blood of vessels D, µm L, µm viscosity µ, cP 1 1 108.0 726.9 3.7 2 2 84.7 549.6 3 4 66.4 4 8 5 Level Number Diameter Length Blood of vessels D, µm L, µm viscosity, µ, cP 1 1 147.0 1036.2 4.2 3.6 2 2 115.3 783.4 4.0 415.5 3.4 3 4 90.4 592.3 3.9 52.1 314.1 3.2 4 8 70.9 447.8 3.7 16 40.8 237.5 2.9 5 16 55.6 338.5 3.5 6 32 32.0 179.5 2.7 6 32 43.6 255.9 3.2 7 64 25.1 135.5 2.4 7 64 34.2 193.5 3.0 8 128 19.7 102.6 2.1 8 128 26.8 146.3 2.7 9 256 15.4 77.6 1.8 9 256 21.0 110.6 2.3 10 512 12.1 58.7 1.5 10 512 16.5 83.6 2.0 11 1024 9.5 44.3 1.2 11 1024 12.9 63.2 1.7 12 2048 7.4 33.5 0.9 12 2048 10.1 47.8 1.4 13 4096 5.8 25.3 2.5 13 4096 7.9 36.1 1.1 14 8192 5.1 21.7 4.2 14 8192 6.2 27.3 2.8 cap. 32768 5.0 500.0 4.6 The dashed line shows the division of the 29 vascular levels from [59] into the arterioles, capillaries and venules compartments in our model, as depicted in Figure 2.2. retinal blood flow and the control value of the pressures P in , P 1,in , and P out , P 5,out (defined below). Pressure The control value of the input pressure for this model is chosen to be two-thirds of the MAP measured at the level of the brachial artery, where the factor two-thirds 15 Table 2.4. Computed values of intraluminal blood pressure at the control state. Value, Pressure Value, Pressure mmHg mmHg P in 62.22 P 3,4 24.25 P 1,in 46.85 P4 22.13 P1 43.92 P 4,5 20.01 P 1,pre 40.99 P 5,post 19.80 P 1,post 40.85 P 5,pre 19.76 P 1,2 40.19 P5 18.84 P2 36.09 P 5,out 17.92 P 2,3 32.00* P out 14.00** P3 28.13 *This value chosen based on previous studies [45, 71]. **This value is chosen to be less than P 5,out and higher than the jugular venous pressure (' 4 − 6 mmHg [72]). accounts for the drop in pressure between the brachial artery and the eye [73]. Since SP = 120 mmHg and DP = 80 mmHg, it follows that P in = 62.2 mmHg, as reported in Table 2.4. The control pressure between arterioles and capillaries, P 2,3 , is set at 32 mmHg [45, 71]. The control pressures at all the other nodes of the network are computed using Ohm’s law and the control values of the resistances, namely P i = P j − Ri,j Q, where the subscripts i and j indicate any two consecutive nodes and Ri,j indicates the resistance between them. The control value of the outlet pressure is chosen to be P out = 14 mmHg, so that it is less than P 5,out and higher than the jugular venous pressure (normally between 4 and 6 mmHg [72]). All control pressures are summarized in Table 2.4. 16 Table 2.5. Computed values of vascular volume and capacitance at the control state. Segment Volume, mL Capacitance, mL/mmHg CRA 2.12 · 10−4 7.22 · 10−7 Arterioles 2.20 · 10−4 7.53 · 10−7 Venules 6.12 · 10−4 1.67 · 10−5 CRV 3.92 · 10−4 1.07 · 10−5 Capacitance The capacitance of a vessel represents its ability to store fluid volume for a given pressure difference between the inside and outside of the vessel. The capacitance of vascular compartments can be computed as the product of vascular volume (V ol) and distensibility (Dist) [74]. The control values of the CRA and CRV volumes 2 2 Lcra /4 and V olcrv = πDcrv Lcrv /4. Analogously, the are computed as V olcra = πDcra control values of the volumes of arterioles, capillaries, and venules are computed by using the diameters and lengths reported in Table 2.3. The distensibilities of the CRA and the retinal arterioles are assumed to be equal to those of the cerebral arteries [74], namely Distcra = Dist2 = 34.12·10−4 mmHg−1 . The distensibilities of retinal venules and the CRV are assumed to be eight times larger than those of the arteries, namely Distcrv = Dist4 = 27.30 · 10−3 mmHg−1 [74]. The control volumes and capacitances are reported in Table 2.5. 2.1.3 Time-profiles of the input and output pressure waves The time profile of Pin and Pout at the control state are determined through an inverse problem based on color Doppler imaging measurements of blood velocity in the CRA and CRV, as shown in Figure 2.3a. The centerline blood velocity V (t) in the 17 CRA and CRV is given by V (t) = 8Q(t)/πD2 , where Q(t) = ∆P (t)/R(t) and ∆P (t) is the pressure drop along the resistor. Thus, equation (2.5) can be used to determine the time profiles Pin (t) and Pout (t) that give the CRA blood velocity profiles shown in Figure 2.3b when the system is at its control state. The control profiles of P in (t) and P out (t) are shown in Figure 2.3c, and their mean values are P in = 62.2 mmHg and P out = 14 mmHg. At the control state, the maximum and minimum values of P in (t) are P in,max = 92.2 mmHg and P in,min = 39.6 mmHg, and their ratios with respect to SP and DP are P in,max /SP = 0.7683 and P in,min /DP = 0.4945. These ratios are used to scale the Pin (t) profile in the simulations of clinical conditions where SP and DP are different from their control state values. (b) (a) (c) 10 100 90 pressure [mmHg] velocity [cm/s] CRA 5 0 0.2 0.4 0.6 0.8 in 70 60 50 40 30 Pout 20 CRV −5 0 P 80 1 time [s] 10 0 0.2 0.4 0.6 0.8 1 time [s] Figure 2.3. (a) Typical color Doppler imaging (CDI) measurement of the blood velocity in the CRA and CRV. (b) Time profiles of the blood velocity in the CRA and CRV at the control state. (c) Time profile of the inlet and outlet pressures P in (t) and P out (t) at the control state. 2.1.4 IOP-induced stress in the lamina cribrosa The lamina cribrosa is modeled as a nonlinear, homogeneous, isotropic, elastic circular plate of radius rlc and finite thickness hlc , satisfying the equilibrium equation, ∇ · S = 0, where S = λlc (σe )tr(E)I + 2µlc (σe )E (2.7) 18 is the stress tensor, h i T T ∇v + (∇v) + (∇v) ∇v E= (2.8) 2 is the Green-Saint Venant strain tensor, v is the displacement vector, µlc is the shear modulus, λlc = µlc Elc − 2µlc 3µlc − Elc (2.9) is the Lamé’s parameter, and Elc is the Young’s modulus. The elastic parameters λlc and µlc vary with the effective stress σe , as described by Guidoboni et al. [48, 75] and in Table 2.7. Table 2.6. Radius and thickness values used in the mathematical model for the lamina cribrosa and sclera Parameter Value Source Radius of the lamina cribrosa rlc , mm 0.75 [76] Thickness of the lamina cribrosa hlc , mm 0.2 [62, 63] Radius of the sclera rs , mm 12 [62, 77] Thickness of the sclera hs , mm 1 [63, 77] Table 2.7. Model values of Young’s modulus Elc and shear modulus µlc for the lamina cribrosa as a function of the effective stress σe , as in [78, 79]. Elc , MPa µlc , MPa Range of σe , MPa 0.358 0.12 0.000 ≤ σe < 0.008 0.656 0.22 0.008 ≤ σe < 0.0015 1.818 0.61 σe ≥ 0.015 19 The problem is solved in cylindrical coordinates, as shown in Figure 2.4, with the boundary conditions Sn = −IOPn at ζ = hlc /2, Sn = −RLTpn at ζ = −hlc /2, nSn = T and vζ = 0 at s = rlc , where T is the scleral tension computed via Laplace’s law IOP rs , 2hs where rs and hs are the radius and thickness of the sclera, respectively. T = IOP ζ ζ θ IOP T s s hlc (2.10) RLTp rlc RLTp Figure 2.4. Schematic representation of geometry and boundary conditions of the elasticity problem for the lamina cribrosa. The anterior surface (ζ = hlc /2) is subject to the intraocular pressure (IOP), while the posterior surface (ζ = −hlc /2) is subject to the retrolaminar tissue pressure (RLTp). The lateral surface (s = rlc ) is connected to the sclera and experiences the scleral tension T . Geometrical and mechanical properties of the lamina and sclera are summarized in Tables 2.6 and 2.7; the numerical solution of the elastic problem via finite elements is described in Guidoboni et al. [48]. The radial component of the normal stress Sss at s = 0 resulting from the solution of this elastic problem is assumed to be the external pressure, PLC (ζ), acting on the intralaminar segments of the CRA and CRV, which are assumed to pierce the lamina in its center [48]. 2.1.5 Variable passive resistance In this Section, the mechanical response of the vessels to changes in the value of transmural pressure difference is modeled. 20 Consider a straight cylinder Ω with cross-section Σ (with general shape, thus not necessarily cylindrical) and length L in a Cartesian coordinate system (x, y, z), as shown in Figure 2.5; that is Ω = Σ×(0, L). Let ex , ey and ez be the unit vectors in the x, y and z (axial) direction, respectively. Let’s remark that an equivalent formulation could be obtained using other coordinate restrictions (e.g. polar coordinates). x y z z=0 Σ z=L Figure 2.5. Representative cylinder Ω in the Cartesian coordinate system (x, y, z) with cross-section Σ, of length L. Let us assume that the motion inside the cylinder Ω can be described by the stationary Navier-Stokes equations:    ∇ · u = 0, in Ω     1 µ   (u · ∇)u = − ∇p + ∆u, in Ω ρ ρ   u = 0, on ∂Σ × (0, L)       p(z = 0) = p0 , p(z = L) = pL , (2.11) where u is the fluid velocity, p is the fluid pressure, ρ is the fluid density and µ is the fluid dinamic viscosity, ∇· is the divergence operator, ∇ is the gradient operator, ∆ is the Laplacian operator and the boldface variables are vectorial quantities. We remark that ρ and µ are assumed to be given positive constants. The third equation in (2.11) represents the no-slip boundary condition. In particular this boundary condition could also be interpreted as    u · n = 0, on ∂Σ × (0, L) Σ Σ u=   uz = 0, on ∂Σ × (0, L), (2.12) 21 where u = uΣ + uz ez , with uΣ = ux ex + uy ey , and nΣ is the outward normal vector to ∂Σ × (0, L). The value of the pressure is assigned at the inlet and at the outlet of the cylinder in (2.11). Starting from (2.11), we will perform a model reduction using the following assumptions: • body forces and mass sources are absent; • p is constant on each Σ ⇒ p = p(z); • uz = u(z)f (Σ), where u(z) represents the average axial velocity on the crosssection Σ and f (Σ) is an appropriate shape function; • axial motion is predominant, thus |uΣ |  uz . A model reduction is performed by integrating the equations in system (2.11) on the cross-section Σ. Although not applied here, it is worth mentioning that a more sophisticated analysis can be obtained by integrating over the infinitesimal volume Σ × (z − dz, z + dz)Z with dz → 0, as shown in [80]. Define Q(z) = Z uz dσ = u(z)A(z) , and A(z) = Σ 1 dσ to be the volumetric Σ average flow rate and the cross-sectional area, respectively. Note that the first equality Z f (Σ) dσ = A(z). together with the definition of uz implies that Σ We begin considering the first equation in (2.11), and we split the divergence operator is into its transversal and axial components to emphasize the contribution of the different directions of motion ∇ · u = ∇Σ · uΣ + ∂uz , ∂z with ∇Σ · uΣ = ∂ux ∂uy + . ∂x ∂y (2.13) 22 Then, integrating over the cross-section, we obtain: Z Z ∇ · u dσ = 0= Σ Z ∇Σ · uΣ dσ + Σ 1 = Σ Z     u · n dγ +  Σ Σ  Z ∂Σ d dz Z d = dz Z 2 = ∂uz dσ ∂z Z uz dσ − Σ (u · nΣ )uz dγ ∂Σ Z u(z)f (Σ) dσ Σ d = u(z) dz 3 Σ ∂uz dσ ∂z   Z       − (u · n )u dγ − (u e · n )u dγ   Σ Σ z z z Σ z   ∂Σ ∂Σ Z f (Σ) dσ = Σ  dQ(z) d u(z)A(z) = , dz dz (2.14) where, in step 1, the contribution of the cross-sectional component of the operator is moved to the boundary of Σ using the divergence theorem and then it is eliminated due to the boundary conditions; in step 2, the partial derivative with respect to z is moved outside of the integral using Leibniz theorem and it is changed to a total derivative; in step 3, u(z) is moved outside of the integral since the integration does not depend on the axial coordinate z. The advective term of the momentum equation in (2.11) can be written as:  u · ∇u = ∇ · (u ⊗ u) −  u∇ · u,  (2.15) where the operator ⊗ represents the outer product. Integrating the second equation in (2.11) on the cross-section Σ and applying the identity (2.15), we obtain: Z Z Z µ 1 ∇ · (u ⊗ u) dσ + ∇p dσ − ∆u dσ = 0. ρ Σ ρ Σ Σ (2.16) Since the contribution of axial motion is predominant and pressure is assumed to be constant on each cross-section Σ, it is possible to focus only on the z−component of equation (2.16), namely Z Z Z 1 ∂p µ ∇ · (uz ⊗ u) dσ + dσ − ∆uz dσ = 0. ρ Σ ∂z ρ Σ Σ (2.17) 23 The first integral in (2.17) becomes Z Z Z ∇ · (uz ⊗ u) dσ = ∇Σ · (uz uΣ ) dσ + Σ Σ Σ ∂u2z dσ ∂z  ∂ u2 (z)f 2 (Σ) = uz uΣ · nΣ dγ + dσ ∂z ∂Σ Σ Z Z  du2 (z)  2  =  f 2 (Σ) dσ uz u Σ · nΣ dγ + ∂Σ dz Σ Z 2 1 du (z) def , where β = = βA(z) f 2 (Σ) dσ, dz A(z) Σ 1 Z Z (2.18) where in step 1 the first integral on the right hand side is moved to the boundary using the divergence theorem, and in step 2 the boundary condition in (2.11) is applied. The second integral in (2.17) becomes Z Z 1 ∂p A(z) dp 1 dp dσ = dσ = . ρ Σ ∂z ρ dz Σ ρ dz (2.19) The third integral in (2.17) becomes µ − ρ Z Z 2 µ ∂ uz ∆Σ uz dσ − dσ ρ Σ ∂z 2 Σ *' 0  Z Z 2  µ µ ∂ u 1 z  dσ =− ∇Σ · (∇Σ uz ) dσ −  2 ρ Σ ρ ∂z   Σ Z  µ = − u(z) ∇Σ · ∇Σ f (Σ) dσ ρ Σ Z µ = − u(z) nΣ · ∇Σ f (Σ) dγ ρ ∂Σ Z µ def nΣ · ∇Σ f (Σ) dγ, = Kr u(z), where Kr = − ρ ∂Σ µ ∆uz dσ = − ρ Σ Z (2.20) where in step 1 the contribution of diffusion in the axial direction is neglected. Combining (2.14), (2.18), (2.19), and (2.20), the following system of equations is obtained:   dQ(z)   = 0, in (0, L) dz 2 du (z) A(z) dp   + + Kr u(z) = 0, in (0, L).  βA(z) dz ρ dz (2.21) 24 Due to the low Reynolds number (Re ' 4 in the CRA), the second equation of (2.21) becomes 2   A(z) dp (z)  du  βA(z) + + Kr u(z) = 0,   dz ρ dz (2.22) Q(z) 1 dp + 2 = 0. Kr ρ dz A (z) (2.23) or equivalently The value of Kr depends on the specific profile f (Σ) that is chosen a priori for the fluid 8πµ motion. In the case of a parabolic velocity profile (Poiseuille’s flow) Kr = [80,81]. ρ Equation (2.23) involves the volumetric flow Q(z) and the pressure gradient along dp the tube and, consequently, the vascular resistance (Ohm’s Law). However, the dz problem of how to introduce the contribution of the transmural pressure difference to the resistance is not well-defined. The problem is closed by introducing the tube law [81–84], P (α) = p(z) − pe . Kp (2.24) A(z) represents the ratio between the cross-sectional area Aref A(z) and the reference cross-sectional area Aref , p(z) represents the pressure inside the def In equation (2.24), α = tube, and pe represents the external pressure acting on the tube, respectively. Thus, equation (2.24) relates the effect of transmural pressure difference (i.e. p(z) − pe ) with the mechanical response of the vessel. The value of the constant Kp is chosen in agreement with the value in [81, 84] for a linear elastic tube. E def Kp = 12(1 − ν 2 ) h rref !3 , (2.25) where E, ν and h are the tube Young’s modulus, Poisson’s ratio and wall thickness respectively. Notice that the reference state corresponds to p − pe = 0 and a circular cross-section, 2 and thus Aref = πrref where rref is the reference state radius of the tube. If p−pe > 0, the cross-section will dilate, corresponding to α > 1, and if p−pe < 0, the cross-section will constrict, corresponding to α < 1. 25 The explicit form of P (α) depends on the physiological and mechanical properties of the considered tube. Arterial walls are thicker than venous walls, and the pressure inside of an arterial branch is usually higher than in the corresponding venous branch. Therefore arteries can be described as compressible tubes that will deform but will unlikely collapse. Veins can be described as collapsible tubes which may be subject to a negative transmural pressure difference and are more likely to collapse. The phenomenon of venous collapse is known as Starling resistor. rref + η rref pe p Σ Figure 2.6. Radial displacement on cross-section Σ in the case of a compressible tube for a positive transmural pressure difference p − pe > 0. For the case of a compressible tube, the cross-section is assumed to remain circular and the radial displacement η is computed using Laplace’s law 2 (p − pe ) (1 − ν 2 )rref . η= Eh (2.26) Some algebraic manipulations are then necessary to introduce P (α) in the equation. 2 (1 − ν 2 )rref (p − pe ) h2 12(1 − ν 2 ) η= = Eh 12rref E  rref h 3 (p − pe ) h2 p − pe h2 p − pe h2 rref = = rref = P (α)rref = P (α), 2 2 12rref Kp 12rref Kp 12rref kL  2 rref where kL = 12 . h η Thus P (α) = kL . In addition rref A = πr2 = π(rref + η)2 , p rref + η = A/π, 2 Aref = πrref , p rref = Aref /π, (2.27) (2.28) 26 and η rref p A/π rref + η rref − 1 = α1/2 − 1. = − =p rref rref Aref /π (2.29) So P (α) = kL (α1/2 − 1). p − pe ε−band 1 α Figure 2.7. Graph of tube law experimental results for a collapsible tube, adapted from [81–84]. The cross-section shape changes from circular to elliptic to highly collapsed as the transmural pressure difference p − pe decreases. For the case of a collapsible tube, when α < 1, the corresponding tube law is P (α) = 1 − α−3/2 . This formula is based on experiments [81, 84] and describes the Starling resistor effect, graphically depicted in Figure 2.7. When the transmural pressure difference is positive, the vessel is dilating and thereby maintaining a circular cross-section. However, when the transmural pressure difference is negative, the cross-section becomes elliptical. In this region, known as ε−band, as shown in Figure 2.7, small changes in transmural pressure difference yield large changes in the cross-sectional area (and in the vascular resistance). As the transmural pressure difference becomes more negative, the vessel continues to collapse until two opposite sides of the inner wall touch each other. In this region, the vessel is more resistant and thus large changes in the transmural pressure difference yield small changes in the cross-sectional area (and in the vascular resistance). 27 The modeling choices adopted for P (α) for a compressible and collapsible tube are summarized below. Compliant tube P (α) = kL (α1/2 − 1),    1 − α−3/2 , α≤1 Collapsible tube P (α) =   kL (α1/2 − 1), α > 1. The last step necessary to close the problem is to use the tube law to find an explicit formulation for the vascular resistance. The case for the collapsible tube will be shown in detail. The computations for compressible tubes are the same as those for the collapsible tubes when α > 1. First, the tube law can be solved for α:    " #−2/3       p(z) − p     e −2/3  ,   , α≤1   1 − P (α)  1− Kp α= = " #2       2   P (α) p(z) − p   e     α>1 +1 ,  k +1 ,   Kp kL L α≤1 (2.30) α > 1. Then, equation (2.23) is integrated from z = 0 to z = L (note that the first equation in (2.21) implies that Q does not depend on z) to obtain an equation for the resistance: Q 1 dp + 2 =0 Kr ρ Zdz A (z) Z L L 1 dp 1 dz = −Q dz 2 Kr ρ 0 dz A (z) 0 Z L pL − p0 Q =− 2 α−2 dz Kr ρ Aref 0 Z p0 − pL Kr ρ L −2 = 2 α dz Q A Z Lref 0 Kr ρ R= 2 α−2 dz. Aref 0 (2.31) 28 Finally the expression for α (2.30) is substituted into this expression for the resistance to obtain        R=       Kr ρ A2ref Z Kr ρ A2ref Z L " 1− 0 0 L " p(z) − pext Kp p(z) − pext Kp kL #4/3 " #4/3 p̂ − pext Kr ρL 1− , α≤1 2 Aref Kp #−4 " #−4 Kr ρL p̂ − pext +1 +1 dz ' 2 , α > 1, Aref Kp kL (2.32) dz ' where the value of the pressure p(z) is approximated by the average pressure p̂ inside p0 + p L the vessel, namely p ' p̂ = . 2 The compliant tube case gives the analogous result " #−4 Kr ρL p − pext +1 R= 2 . Aref Kp kL (2.33) It is also possible to check that the modeling approach used in the compliant tube case is equivalent to Poiseuille’s Law " #−4  −4 8πµρL P (α) Kr ρL p − pext +1 = 2 4 +1 R = 2 Aref Kp kL ρπ rref kL " #−4 " #−4 8µL η 8µL η + r (2.27) ref = +1 = 4 = 4 πrref rref πrref rref " #−4 8µL r 8µL 128µL = 4 = = , 4 πrref rref πr πD (2.34) where r and D represent the tube radius and diameter, respectively. The compliant tube case is used to compute the resistances R1c and R1d (Figure 2.1), while the collapsible tube case formula is used to compute the resistances R4a , R4b , R5a and R5b . The value of the external pressure pe varies depending on the region in which the resistance is positioned (pe = IOP in the intraocular region and pe = PLC (ζ) in the translaminar region). The physical parameters of the segments in the venules compartment used in the computation of R4a and R4b are chosen as follows: the values of the wall-to lumen-ratio at control state and the Poisson’s ratio are h/D = 0.05 and ν = 0.49, respectively, 29 where h represents the thickness of the vessel wall and D is the diameter. These values are the same as the values in the CRV. The value of the Young’s modulus is E = 0.066 MPa. 2.1.6 Active variable resistance The resistances R2a and R2b are modeled through a phenomenological description of blood flow autoregulation, following the method utilized by Lakin et al. in the context of cerebral blood flow [74]. The absence of autoregulation is modeled by keeping the arteriolar resistances R2a and R2b constantly equal to their control value R2a and R2b , respectively, despite changes in OPP, while functional autoregulation is modeled by allowing the arteriolar resistances to vary according to the formula:     cL + cU exp K QnoAR (OPP) − Q − ĉ    , (2.35) R2a = R2b = R2a  1 + exp K QnoAR (OPP) − Q − ĉ where cL = 2.5 · 10−3 and cU = 1.18 determine the lower and upper bounds for the resistance, K = 3.45 · 104 s/mL determines the sensitivity of vessel resistance to changes in OPP, QnoAR is the total retinal blood flow predicted by the model in the absence of autoregulation for a given OPP and ĉ = ln(cU −1)−ln(1−cL ) ensures that R2a = R2a and R2b = R2b at the control state. The numerical values of cL , cU and K have been chosen to yield an increase in resistance over the range of OPP depicted in Figure 2.8a. Model simulations of normalized flow in the case of functional and absent autoregulation are reported in Figure 2.8b for various OPP values (attained by setting IOP = 15 mmHg and varying MAP between 63 and 123 mmHg). Figure 2.10c shows that the model predicted values of retinal blood flow in the case of functional autoregulation are consistent with clinical measurements for various MAP levels; this suggests that formula (2.35) and its related assumptions are appropriate modeling choices for retinal blood flow autoregulation. 30 (a) (b) 8000 7000 1.3 R̄2a 1.2 normalized flow R2a [mmHg s/ml] 6000 5000 4000 3000 2000 1.1 1 0.9 0.8 functional autoregulation 1000 0.7 absent autoregulation 0 30 40 50 OPP [mmHg] 60 70 0.6 30 40 50 60 70 OPP[mmHg] Figure 2.8. Arteriolar resistance (a) and normalized flow (b) for model simulations with functional (solid line) or absent (dashed line) autoregulation for various ocular perfusion pressures (OPP) attained by setting IOP = 15 mmHg and varying MAP. 2.2 Solution algorithm In this Section we illustrate the solution procedure for the nonlinear problem in system (2.5). The nonlinearities of the system are due to the variable passive resistances that are highly nonlinear functions of the nodal pressures in Figure 2.1. To address these numerical difficulties, the problem is linearized and solved with an iterative algorithm using the fixed point method as described below: Given the values of IOP, RLTp, SP, DP, and an initial guess for each nodal pressure in Figure 2.1 Pj0,0 , def  for j ∈ H = (in), (1, in), (1), . . . , (5), (5, out), (out) . (2.36) Then, let n = 0, . . . , N be the nodes of the time discretization of a cardiac cycle: • Step 1: 31 (intraocular ζ portion) hLC /2 IOP lamina central retinal vessels PLC (ζ) cribrosa IOP PLC (ζ) cribrosa (translaminar portion) −hLC /2 RLTp lamina RLTp (retrobulbar portion) Figure 2.9. Schematic representation of the translaminar portion of the central retinal vessels. n to be consistent with the given values of SP and DP – Update Pin – Solve thelamina cribrosa problem in [48] to obtain the value of PLC (ζ),  hLC hLC with ζ ∈ − , 2 2 – Set R2a    R , if AR is absent 2a = R2b =   (2.35), if AR is functional; (2.37) for k ≥ 0 • Step 2: For n = 1, . . . , N n,k – Compute R1c using a modified version of formula (2.33) that includes the variation along the coordinate ζ of the value of PLC as shown in Figure 2.9. #−4 Z hLC /2 " K ρ P̂ − P (ζ) r LC n,k +1 dζ, (2.38) = 2 R1c Aref −hLC /2 Kp kL where P̂ =  n−1,k P1,pre − n−1,k P1,post  and the values of Kr , Aref , Kp and kL are computed for the translaminar portion of the central retinal artery 32 n,k using formula (2.33) – Compute R1d n,k R1d Kr ρ = 2 Aref " P̂ − IOP +1 Kp kL #−4 , (2.39)   n−1,k n−1,k where P̂ = P1,post − P1,2 and the values of Kr , L, Aref , Kp and kL are computed for the intraocular portion of the central retinal artery n,k n,k – Compute R4a and R4b using a modified version of formula (2.32) to ac- count for the architecture from Takahashi et al. [59]  " #4/3   Kr ρL P̂ − IOP   1− ,   A2 Kp ref n,k n,k #−4 " R4a = R4b = 0.5 13  X  Kr,i ρLi P̂i − IOP   +1 ,   A2ref,i Kp,i kL,i i=1 if P̂ − IOP ≤ 0 if P̂ − IOP > 0, (2.40) where P̂ = P4n−1,k and the values of Kr , L, Aref , Kp and kL are computed for an equivalent venular compartment obtained by lumping the venular levels from 1 to 13 [59], while the values of P̂i , Kr,i , Li , Aref,i , Kp,i and kL,i are computed for each one of the 13 venular compartments in [59]. This approach was adopted to eliminate the oscillations present in the system due the collapsibility of the 13 levels of venules from [59] n,k – Compute R5a        n,k R5a =       where P̂ =  using formula (2.32) #4/3 " Kr ρL P̂ − IOP 1− , 2 Aref Kp " #−4 Kr ρL P̂ − IOP +1 , A2ref Kp kL n−1,k P4,5 − n−1,k P5,post  if P̂ − IOP ≤ 0 (2.41) if P̂ − IOP > 0, and the values of Kr , L, Aref , Kp and kL are computed for the intraocular portion of the central retinal vein 33 n,k using a modified version of formula (2.32) that includes the – Compute R5b variation along the coordinate ζ of the value of PLC as shown in Figure 2.9  " #4/3   K ρ P̂ − P (ζ)  r LC  1− , if P̂ − PLC (ζ) ≤ 0   A2 Kp ref " #−4 Υ(ζ) = (2.42)   P̂ − P (ζ) K ρL LC r   +1 , if P̂ − PLC (ζ) > 0   A2 Kp kL ref and n,k R5b Z −hLC /2 Υ(ζ) dζ, = (2.43) hLC /2 where P̂ =  n−1,k P5,post − n−1,k P5,pre  and the values of Kr , Aref , Kp and kL are computed for the translaminar portion of the central retinal vein; • Step 3: Solve the linearized version of system (2.5) for Pjn,k , j = 1, 2, 4, 5    n  d(P1n,k − RLTp) Pin − P1n,k P1n,k − P2n,k    C = − 1  n,k n,k  dt Rin + R1a R1b + R1c + R1d + R2a      n,k n,k n,k  P1 − P2 P2n,k − P4n,k d(P2 − IOP)    = C −  2 n,k n,k n,k dt R1b + R1c + R1d + R2a R2b + R3a + R3b + R4a   P2n,k − P4n,k P4n,k − P5n,k d(P4n,k − IOP)    = − C4  n,k n,k n,k n,k  dt R2b + R3a + R3b + R4a R4b + R5a + R5b + R5c      n,k n,k n,k n,k n  P4 − P5 P − Pout d(P5 − RLTp)    = n,k − 5 ; C  n,k n,k  5 dt R5d + Rout R4b + R5a + R5b + R5c (2.44) • Step 4: Compute the value of the retinal blood flow Qn,k and update the value of the pressures Pjn,k j ∈ H using Ohm’s law; • Step 5: If k ≥ 1 check for the periodicity of the solution of system (2.44) ! maxn=0,...,N |Pjn,k − Pjn,k−1 | < φ, (2.45) max j=1,2,4,5 maxn=0,...,N |Pjn,k | where φ is a given tolerance; • Step 6: If (2.45) is satisfied, then Pjn,k j ∈ H, n = 0, . . . , N is the pressure solution at the nodes of the system in Figure 2.1 for a complete cardiac cycle. If (2.45) is not satisfied then set Pj0,k+1 = PjN,k j ∈ H and return to Step 2. 34 2.3 Results Six different cases corresponding to different clinical conditions are simulated us- ing the model (2.5). The six cases represent patients with high, normal and low arterial blood pressure (HBP-, NBP-, LBP-), with functional or absent blood flow autoregulation (-wAR, -woAR), as summarized in Table 2.8. Systolic/diastolic arterial blood pressures for cases HBP-, NBP- and LBP- are assumed to be 140/90 mmHg, 120/80 mmHg and 100/70 mmHg, respectively. Functional autoregulation is simulated by allowing arteriolar resistances R2a and R2b to vary according to equation (2.35). To simulate the absence of autoregulation, arteriolar resistances R2a and R2b do not change and are set equal to their control values R2a . and R2b . Table 2.8. Representative cases for model simulations. The cases represent individuals with normal, high and low arterial blood pressure (NBP-, HBP-, LBP-) and with functional or absent blood flow autoregulation (-wAR, -woAR). Case Acronym Description High BP HBP- Systolic/diastolic blood pressure = 140/90 mmHg Normal BP NBP- Systolic/diastolic blood pressure = 120/80 mmHg Low BP LBP- Systolic/diastolic blood pressure = 100/70 mmHg Functional AR -wAR Variable arteriolar resistances, see equation (2.35) Absent AR -woAR Arteriolar resistances held fixed at control values 2.3.1 Model validation The model predicted mean values of blood velocity and flow along the retinal vascular network obtained for the NBPwAR case (which represents a healthy clinical condition) are compared in Figures 2.10a and 2.10b with measurements obtained by Garcia et al. [85] and Riva et al. [86] using bidirectional laser Doppler velocimetry 35 in healthy individuals. The pressures Pi at each node of the network are computed from (2.5); the velocity and flow in arterioles, capillaries and venules are determined from these pressures and the calculated resistances R. The model predicted values of total retinal blood flow for IOP = 15 mmHg and MAP between 65 and 115 mmHg are compared with clinical data [49, 55–57] in Figure 2.10c. The model predicted values are normalized with respect to the control state, and the measured data are normalized with respect to their reported baseline value. (a) (b) 5 1.5 Arterial − Venous Arterial − Venous 1.4 20 3.5 3 2.5 2 1.5 1 normalized flow 4 flow [µl /min] velocity [cm/s] 4.5 (c) 25 15 10 5 0.5 0 200 1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 100 0 100 vessel diameter [ µm] 200 0 200 100 0 100 vessel diameter [ µm] 200 0.5 70 80 90 100 110 MAP [mmHg] Figure 2.10. (a) Comparison of model predicted values (solid line) of blood velocity with measured data (open triangles [85], black dots [86]) for vessels of various diameter. (b) Comparison of model predicted values (solid line) of volumetric blood flow with measured data (open triangles [85], black dots [86]) for vessels of various diameter. (c) Comparison of model values (solid line) of normalized total retinal blood flow with measured data (black dots [55], open triangles [57], open squares [57], open diamonds [56], plusses [49]) for various mean arterial pressures (MAP). In all model simulations, IOP = 15 mmHg. 2.3.2 Theoretical investigations The mathematical model is used to investigate the effects of blood pressure and blood flow autoregulation on the IOP-induced hemodynamic changes in total retinal blood flow, CRA blood velocity and intraluminal blood pressure along the retinal vasculature. The model predictions for total retinal blood flow (computed as the 36 time-average of Q(t) over a cardiac cycle), peak systolic velocity (PSV, maximum value of centerline blood velocity along a cardiac cycle), end diastolic velocity (EDV, minimum value of centerline blood velocity along a cardiac cycle) and resistive index (RI = (PSV-EDV)/PSV) in the CRA are compared for cases NBPwAR, NBPwoAR, HBPwAR, HBPwoAR, LBPwAR and LBPwoAR in Figure 2.11. functional autoregulation retinal blood flow [µl/min] (a) CRA − EDV [cm/s] (c) 50 3 10 30 0.9 2 20 HBPwAR 5 0.8 1 10 NBPwAR LBPwAR 20 30 40 0 20 30 40 0 20 30 40 IOP [mmHg] IOP [mmHg] IOP [mmHg] (e) (f) (g) 0.7 20 30 40 IOP [mmHg] (h) 50 1 40 3 10 30 0.9 2 20 HBPwoAR 5 0.8 1 10 0 CRA − RI (d) 1 40 0 absent autoregulation CRA − PSV [cm/s] (b) NBPwoAR LBPwoAR 20 30 40 IOP [mmHg] 0 20 30 40 IOP [mmHg] 0 20 30 40 IOP [mmHg] 0.7 20 30 40 IOP [mmHg] Figure 2.11. Model predicted values of total retinal blood flow; peak systolic velocity in the central retinal artery (CRA-PSV); end diastolic velocity in the CRA (CRA-EDV); and resistivity index in the CRA (CRARI) as IOP varies between 15 and 45 mmHg for theoretical patients with low, normal or high blood pressure (LBP-, NBP-, HBP-) and functional or absent blood flow autoregulation (-wAR, -woAR). Figure 2.12 reports the IOP ranges in 15 − 30 mmHg for which the values of total retinal blood flow (Q) and the PSV, EDV and RI in the CRA are within ±3% of their control values, as predicted by the mathematical model for six theoretical patients. In the low blood pressure cases (LBPwAR and LBwoAR), the values of Q and EDV are never within 3% of their control values. The small (3%) change in blood flow is chosen as a rough threshold for quantifying the regime when autoregulation is successful. 37 IOP [mmHg] 30 25 20 LBPwoAR LBPwAR NBPwoAR NBPwAR HBPwoAR PSV EDV RI PSV EDV RI Q PSV EDV RI Q PSV EDV RI Q PSV EDV RI Q PSV EDV RI Q Q 15 HBPwAR Figure 2.12. Model predicted ranges of IOP in 15−30 mmHg for which the values of total retinal blood flow (Q, black), PSV (blue), EDV (red) and RI (green) are within ±3% of their control state values for six theoretical patients with low, normal or high blood pressure (LBP-, NBP-, HBP-) and functional or absent blood flow autoregulation (-wAR, -woAR). Total retinal blood flow As shown in Figures 2.11 and 2.12, the model predicted average blood flow over a cardiac cycle remains relatively constant when IOP= 15 − 25 mmHg in NBPwAR individuals. In the HBPwAR case, relatively constant flow is predicted when IOP= 25 − 30 mmHg. In the remaining cases, the model predicts a slight or nonexistent autoregulation plateau when IOP= 15 − 30 mmHg. CRA blood velocity Although a monotone decrease in PSV and EDV with IOP elevation is predicted in most of the theoretical patients, Figures 2.11b and 2.11c show that the NBPwAR and HBPwAR cases exhibit a slight PSV increase for IOP within the regulating range, namely between 15 and 25 mmHg for NBPwAR and between 25 and 30 mmHg for 38 HBPwAR. For all six cases, the model predicts an increase in the CRA resistivity index with IOP, as shown in Figures 2.11d and 2.11h. As shown in Figure 2.12, for all the cases except HBP-, the range of IOP for which PSV is within 3% of its control values is larger than the range for EDV. In the HBPwAR, case the PSV value is greater than 3% of its baseline value for all of the simulated values of IOP. Intraluminal blood pressure The model predicted values of blood pressure in the five model compartments of the retinal vasculature for IOP equal to 15, 30 and 45 mmHg are compared in Figure 2.13 for the six representative clinical cases described in Table 2.8. Pronounced changes in pressure are visible in the CRA, arterioles, capillaries and venules with IOP elevation, while the pressure in the CRV remains relatively constant. As expected, when IOP is sufficiently high, the venules collapse. The degree to which the vessels collapse depends on the MAP of the individual. In particular, the model predicts that for an IOP = 45 mmHg, the venules collapse to a lesser extent in HBP- patients than in NBP- and LBP- patients. 2.3.3 Comparison between model predictions and clinical data The mathematical model is used to interpret clinically measured hemodynamic responses to surgical IOP reduction in glaucoma patients and induced IOP elevation in healthy individuals. Surgical IOP reduction on glaucoma patients Galassi et al. [10] and Trible et al. [15] report an increase in CRA blood velocity and a decrease in RI following trabeculectomy (group (a) in Galassi et al. and Trible et al.) and deep sclerectomy (group (b) in Galassi et al.), as indicated by the changes (∆) in mean flow velocity (MFV), PSV, EDV and RI summarized in Table 2.9. Since IOP = 15 mmHg IOP = 30 mmHg IOP = 45 mmHg pressure [mmHg] pressure [mmHg] pressure [mmHg] 39 CRA (P1) arterioles (P 2) capillaries (P3) venules (P4) CRV (P5) 100 100 100 100 100 80 80 80 80 80 60 60 60 60 60 40 40 40 40 20 20 0 0.5 1 20 0 0.5 1 40 20 0 0.5 1 20 0 0.5 1 0 0.5 1 100 100 100 100 100 80 80 80 80 80 HBPwoAR 60 60 60 60 60 NBPwAR 40 40 40 40 40 20 20 20 20 20 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 100 100 100 100 80 80 80 80 80 60 60 60 60 60 40 40 40 40 40 20 20 20 20 0.5 time [s] 1 0 0.5 time [s] 1 0 0.5 time [s] 1 NBPwoAR LBPwAR LBPwoAR 0 100 0 HBPwAR 0.5 1 0.5 1 20 0 0.5 time [s] 1 0 time [s] Figure 2.13. Model predicted values of blood pressure in CRA (P1 ), arterioles (P2 ), capillaries (P3 ), venules (P4 ), and CRV (P5 ) for IOP = 15, 30, and 45 mmHg for six theoretical patients with high (blue), normal (black) and low (green) blood pressure (HBP-, NBP-, LBP-), and functional (solid) or absent (dashed) blood flow autoregulation (-wAR, -woAR). the study by Trible et al. [15] reports MAP values of 107 ± 14.7 mmHg, we use HBP-model predictions that correspond to MAP = 106.7 mmHg as a comparison with this set of clinical data. Also, to be consistent with the data, model simulations assume that IOP = 30 mmHg before surgery and IOP = 15 mmHg after surgery. The model assumptions have been indicated with shaded boxes in Table 2.9. Model predicted values of ∆MFV, ∆PSV, ∆EDV and ∆RI for the HBPwoAR case and the trabeculectomy data reported by Trible et al. [15] differ only by 0.35 cm/s, 0.54 cm/s, 0.56 cm/s and 0.02, respectively. The model predicted values of ∆EDV and ∆RI for the HBPwAR and HBPwoAR cases and the deep sclerectomy data reported by Galassi et al. [10] (group (b)) also show the same qualitative trends. 40 Table 2.9. Changes in peak systolic velocity (PSV), end diastolic velocity (EDV), resistivity index (RI) and mean flow velocity (MFV) in the central retinal artery as IOP is reduced by trabeculectomy (group (a) in Galassi et al. [10] and Trible et al. [15] ), deep sclerectomy (group (b) in Galassi et al. [10]) and model simulations. Model assumptions are reported in shaded boxes. Parameter Study Clinical Studies Mathematical Model [10] Group (a) [10] Group (b) [15] HBPwAR HBPwoAR Pre-op IOP, mmHg 25.27 ± 6.24 24.05 ± 4.27 27.5 ± 8.8 30 30 Post-op IOP, mmHg 9.86 ± 2.10 10.79 ± 2.27 12.5 ± 10.5 15 15 Post-op MAP, mmHg − − 107 ± 14.7 106.7 106.7 Post-op MAP, mmHg − − 102.2 ± 11.5 106.7 106.7 ∆MFV, cm/s − − 1.4 1.05 1.75 ∆PSV, cm/s − − 1.2 0.15 1.73 ∆EDV, cm/s 0.40 0.47 0.96 1.32 1.52 −0.06 −0.05 −0.12 −0.11 −0.10 ∆RI Model assumptions are reported in boldface. Dashes indicate that the quantity was not reported in the article. Induced IOP elevation on healthy individuals The model predicted values of PSV and EDV in the CRA for IOP between 15 and 45 mmHg are compared with clinical data [11] in Figure 2.14a. In the study by Harris et al. [11], PSV and EDV measurements were performed at approximately 3 mm behind the optic disc surface, corresponding to the retrobulbar CRA segment in our model. Since the study involved healthy individuals, the clinical data are compared with the NBP-model predictions. Figure 2.14a shows that the PSV values measured by Harris et al. and predicted by our model in the NBP- case fall precisely in the same range. The model-predicted values of the percent change in the CRA mean flow velocity with IOP elevation are compared with the clinical data by Findl et al. [9] in Fig- 41 ure 2.14b. In the study by Findl et al., IOP was raised 20 mmHg above baseline in two steps and then was reduced by 10 mmHg before releasing the suction cup. Measurements were performed at approximately 3 mm behind the optic disc surface, corresponding to the retrobulbar CRA segment in our model. Findl et al. considered healthy individuals and therefore the clinical data are compared with the NBP- model predictions. (a) (b) 0 % change of CRA mean flow velocity CRA blood velocity [cm/s] 12 PSV 10 8 Data − Harris et al (1996) 6 Model − NBPwoAR Model − NBPwAR 4 EDV 2 0 15 20 25 30 35 40 45 IOP [mmHg] −5 −10 −15 −20 −25 −30 −35 −40 Data − Findl et al (1997) Model − NBPwAR −45 Model − NBPwoAR −50 10 20 10 Increase in IOP [mmHg] Figure 2.14. (a) Comparison between peak systolic velocity (PSV) and end diastolic velocity (EDV) measured by Harris et al. [11] in the central retinal artery (CRA) and the model predicted values in the NBPwAR and NBPwoAR cases. (b) Comparison between the percent change in CRA mean flow velocity measured by Findl et al. [9] and the model predicted values in the NBPwAR and NBPwoAR cases. 2.3.4 Trabeculectomy simulations Alterations in retinal hemodynamics following trabeculectomy (TB) have been reported by several clinical studies [10,14,15], although no changes in blood velocity in the CRA where observed after trabeculectomy in other studies [18,25]. This suggests 42 that individual factors might mediate the hemodynamic outcomes of the procedure, including arterial blood pressure and autoregulation [37, 87]. The model results are compared with the clinical data collected in a study involving 30 patients (age 62.7 ± 9.7) who underwent trabeculectomy [88]. In this study, IOP and PSV and EDV in the CRA were measured before and after trabeculectomy in both the healthy eye and the eye on which surgery was performed. Arterial blood pressure was also measured before and after trabeculectomy and did not show significant changes. Table 2.10. Representative cases for trabeculectomy model simulations. The cases represent individuals with normal and high blood pressure (NT-, HT-) and with functional or absent blood flow autoregulation (-wAR, -woAR). The values of pressure considered in the case HT are the same average values reported in [88]. Case Acronym Description High BP HT- Systolic/diastolic blood pressure = 143/84 mmHg Normal BP NT- Systolic/diastolic blood pressure = 120/80 mmHg Functional AR -wAR Variable arteriolar resistances, see equation (2.35) Absent AR -woAR Arteriolar resistances held fixed at control values The values of PSV, EDV, Q in the CRA are computed with the mathematical model for (i) IOP= 28.7 mmHg and 13.1 mmHg, corresponding to IOP levels before and after trabeculectomy; (ii) systolic and diastolic blood pressures SP/DP= 120/80 mmHg and 143/84 mmHg, representing normotensive (NT-) and hypertensive (HT-) subjects; and functional autoregulation (-wAR) or absent autoregulation (-woAR), as reported in Table 2.10. Statistical analysis of the clinical data shows that PSV and EDV increase from 7.81 ± 2.51 cm/s and 2.39 ± 0.94 cm/s to 8.90 ± 2.66 cm/s and 2.77 ± 1.14 after TB, as IOP decreases from 28.7 ± 7.9 mmHg to 13.1 ± 5.8 mmHg (reported in Table 2.11). 43 The mean values of SP/DP in [88] are 143/84 mmHg, as for the HT case of the mathematical model. Table 2.11. Changes in IOP and CRA hemodynamics following trabeculectomy [88]. Pre-Trabeculectomy Post-Trabeculectomy Parameter Healthy Eye Operated Eye Healthy Eye Operated Eye IOP, mmHg 17.9 ± 5.7 28.7 ± 7.9 16.4 ± 5.2* 13.1 ± 5.8* PSV, cm/s 8.73 ± 1.49 7.81 ± 2.51 8.59 ± 2.28 8.90 ± 2.66* EDV, cm/s 2.89 ± 1.28 2.39 ± 0.94 2.75 ± 1.04 2.77 ± 1.14* *Paired Samples t test. Significance level p< 0.05. Additional information on the velocity variations in each patient after trabeculectomy is provided in the scatter plots reported in Figure 2.15, where changes in PSV and EDV are reported as a function of MAP for each patient. While the majority of patients showed an increase in PSV and EDV after trabeculectomy, in agreement with the trend of the average values in Table 2.11, some patients showed a decrease in either PSV or EDV following trabeculectomy. The model predictions for PSV and EDV values in the CRA over a cardiac cycle before and after trabeculectomy are displayed in Figure 2.16. The model predicts an increment in both PSV and EDV for the cases HTwAR, HTwoAR and NTwoAR as IOP is reduced from 28.7 mmHg to 13.1 mmHg. However, for the case NTwAR, the model predicts an increase of EDV and a decrease of PSV values. In particular, Figure 2.16c shows that for almost half of the cardiac cycle the velocity profile before TB is below the value predicted after TB since the difference between the values before and after TB is negative. In Figure 2.17, clinical observations are simulated using the model by calculating the mean volumetric retinal blood flow over a cardiac cycle. The mean blood flow 44 (a) PSV post-TB − PSV pre-TB (b) EDV post-TB − EDV pre-TB Figure 2.15. Scatter plot of the velocity difference in the CRA for the patients in [88]. PSV (a) and EDV (b). (b) IOP=13.1 mmHg 12 CRA velocity [cm/s] CRA velocity [cm/s] 12 10 8 6 4 2 0 0 (c) V post−TB − V pre−TB velocity difference [cm/s] (a) IOP=28.7 mmHg 0.5 time [s] 1 10 8 6 4 2 0 0 0.5 time [s] 1 2.5 2 1.5 PSV 1 EDV 0.5 0 −0.5 0 0.5 time [s] 1 Data − Diliene et al.(2013) Model − NTwoAR Model − NTwAR Model − HTwoAR Model − HTwAR Figure 2.16. Model predictions of: CRA velocity before TB (a); CRA velocity after TB (b). Model predictions of the difference between after and before TB of: CRA velocity, and clinical data (c) [88]. 45 increases as IOP is decreased following trabeculectomy in the HTwoAR and NTwoAR cases, but remains almost constant in the NTwAR case. The HTwAR case shows an intermediate trend. The dotted line in Figure 2.17 represents the IOP value at which we have the intersection between the average retinal blood flow values for the cases NTwAR and NTwoAR. −4 9 x 10 NTwoAR NTwAR 8.5 retinal blood flow [ml/s] HTwoAR HTwAR 8 7.5 7 6.5 6 5.5 5 4.5 12 14 16 18 20 22 24 26 28 30 IOP [mmHg] Figure 2.17. Model predictions of average retinal blood flow in a cardiac cycle for the four theoretical patients. 2.4 Discussion The combination of multiple factors, including IOP, BP and autoregulation, plays a primary role in determining retinal and retrobulbar hemodynamics. It is widely recognized that IOP elevation poses a serious challenge to tissue perfusion, and several animal and human studies, in addition to the results of the theoretical model presented here, have suggested that arterial blood pressure [3, 37–41], and blood flow autoregulation [3, 42–44], are also important factors influencing blood flow. However, 46 the difficulty of isolating and assessing the contributions of arterial blood pressure and blood flow autoregulation in vivo limits the current understanding of their effects on tissue perfusion as IOP varies [37]. This work introduces a mathematical model that can isolate and assess the contributions of arterial blood pressure, blood flow autoregulation and IOP to retinal hemodynamics. 2.4.1 Model validation Every mathematical model is based on simplifying assumptions whose validity needs to be verified by comparing model predictions with data from independent experimental and clinical studies. Figure 2.10 shows that the model predicted values of velocity and flow are consistent with clinical measurements [49,55–57,85,86]. These measurements were used for comparison purposes only and were not used to estimate any of the model parameters listed in Tables 2.1, 2.2, 2.4, 2.5, 2.6, and 2.7. It is important to note that many of the model parameters vary among individuals. For example, blood pressure and elasticity of the vessel walls differ depending on age [89–92], ethnicity [93, 94] and disease [95, 96]; many studies have shown that the geometric and mechanical properties of the lamina cribrosa and sclera undergo changes with age [97, 98], ethnicity [75, 99], and disease [100, 101], and these changes strongly affect the biomechanical response of the optic nerve head tissues to IOP alterations [102–108]. The clinical studies [49,55–57,85,86] used for model comparison did not report the parameter values for the geometric and mechanical properties of the lamina cribrosa and sclera, and therefore the model simulations were performed using the literature-based parameter choices reported in Tables 2.6 and 2.7. Future data on the geometric and structural properties of the retinal vessels collected from populations including a wide age range will help improve upon the current model predictions. Despite the simplified vascular architecture and the literature-based parameter choices, the model predictions of velocity and flow reported in Figure 2.10 are con- 47 sistent with six independent clinical studies, providing evidence that our modeling choices are appropriate and physiologically reasonable. 2.4.2 Theoretical investigations The model simulations suggest that the hemodynamic response of the retinal vasculature to IOP variations differs noticeably among individuals with different blood pressure and functionality of autoregulation. Total retinal blood flow The model predicts that the autoregulation plateau would shift towards higher IOP values as MAP increases, as shown in Figure 2.12. This is in agreement with the study by He et al. [39], who induced IOP elevations on Long-Evan rats with low, moderate, and high MAP levels and found that a higher IOP was needed to attenuate ocular blood flow in animals with higher MAP. The shift in the autoregulation plateau was also predicted by the theoretical work of Arciero et al. [45], whose microcirculation model suggested that autoregulation fails to operate over its expected range of arterial pressure if IOP is increased. CRA blood velocity The model predicts that the blood flow velocity in the CRA does not always decrease with IOP elevation. Figures 2.11b and 2.11c show a slight increase in PSV in the range of IOP values for which autoregulation is achieved, namely between 15 and 25 mmHg for NBPwAR and between 25 and 30 mmHg for HBPwAR. These findings have important clinical implications. For example, the model suggests that a 10 mmHg IOP reduction in an individual with normal blood pressure and functional autoregulation (NBPwAR) would result in a noticeable increase in PSV and EDV in the CRA blood velocity and a noticeable increase in total retinal blood flow if 48 IOP is reduced from 40 to 30 mmHg. However, only minimal hemodynamic changes are predicted if IOP is reduced from 25 to 15 mmHg. These findings could help to explain why some studies [10, 14, 15] report significant hemodynamic changes following trabeculectomy, while others do not [18, 25]. The model also suggests that IOP reductions of approximately 5 mmHg would result in minor hemodynamic changes in all of the cases considered here, which is consistent with clinical data related to topical medications [8,12,13,16,17,20–24,26–28]. However, reductions in IOP greater than 5 mmHg may have a more significant effect on retinal hemodynamics. Intraluminal blood pressure Intraluminal blood pressure is the main driving force of local tissue perfusion and can be measured in retinal vessels by artificially increasing IOP using ophthalmodynamometry [109–112]. However, the level of IOP may have an important impact on the intraluminal pressure in retinal vessels, which is uncovered by the mathematical model. The model predicts that increased IOP induces a significant increase in the intraluminal blood pressure in all vascular compartments upstream of the CRV. This finding is consistent with the experimental observations by Gluksberg and Dunn [113] and Attariwala et al. [114] on live anesthetized cats. A hydraulic feedback mechanism could explain this phenomenon [113, 114]. Veins are more susceptible than arteries to IOP elevation, since veins have thinner walls and lower intraluminal pressure than arteries and, under extreme conditions, act like a Starling resistor and collapse. Thus, as IOP increases, resistance to flow increases in veins more than in arteries, leading to an overall increase in intraluminal pressure upstream of the veins [115, 116]. In addition, the model predicts that IOP elevation would affect intraluminal arterial pressure differently depending on whether blood flow autoregulation is functional or absent. The study by Jonas [110] supports this finding, since correlation coefficients between intraluminal blood pressure measured via ophthalmodynamometry and sys- 49 temic blood pressure were found to be lower in eyes with retinal or orbital diseases than in the control group. 2.4.3 Comparison between model predictions and clinical data A good qualitative and quantitative agreement in the hemodynamic response to IOP changes was found between clinical data and model predictions. This suggests that the mathematical model could be used to anticipate the hemodynamic outcome of clinical IOP modulation on specific patients. Surgical IOP reduction on glaucoma patients The clinical studies by Galassi et al. [10] and Trible et al. [15] suggest that clinically measurable changes in the retinal hemodynamic parameters occur in patients with elevated blood pressure as IOP is reduced from 30 to 15 mmHg. The theoretical analysis in the current study suggests that these changes are less pronounced in HBPwAR patients than in HBPwoAR patients, likely due to the compensatory mechanisms of autoregulation [117]. However, the model also suggests that if a 15 mmHg IOP reduction falls within the regulating range predicted for HBP-patients, the consequent hemodynamic changes would be negligible. This might explain why the trabeculectomy study by Cantor [18] did not observe any significant change in ocular blood flow parameters despite significant IOP reduction. Induced IOP elevation on healthy individuals The ability to autoregulate blood flow has a noticeable effect in individuals with normal blood pressure experiencing induced IOP elevation. Figure 2.14a shows good agreement between clinical data and model predictions in the NBPwoAR case when IOP is increased. Blood velocity measurements were performed immediately after IOP elevation (i.e., before autoregulation took effect), explaining why the NBPwoAR 50 model predictions match the data more closely than the NBPwAR predictions. Consistent with this interpretation, there is better agreement between the data and the NBPwAR case in Figure 2.14b since the blood velocity measurements were performed at least 5 minutes after IOP elevation, leaving time for autoregulation to exert its influence. 2.4.4 Trabeculectomy simulations The clinical data from Table 2.11 show that, on average, trabeculectomy leads to an increase in both PSV and EDV as IOP is reduced. However, Figure 2.15 shows a more complicated scenario in which both increases and decreases in PSV and EDV are possible outcomes of trabeculectomy depending on the patient. The model predictions show how the hemodynamic outcomes following trabeculectomy can vary among patients depending on blood pressure and autoregulatory capacity. For a normotensive patient, the same IOP reduction following trabeculectomy can cause opposite variations in PSV and EDV depending on the status of the autoregulatory mechanisms, as shown in Figure 2.16. For a hypertensive patient the autoregulatory mechanism starts working at a higher level of IOP, which means that when IOP is reduced, the system falls outside of the autoregulatory range. This explains why the change in velocity in the HTwAR case is qualitatively similar to the NTwAR case but quantitatively different in Figure 2.16c. The level of blood pressure and the functionality of AR also influence the model predictions of the retinal blood flow value over a cardiac cycle. When autoregulation is within its working range (NTwAR), flow remains almost constant as IOP is reduced following trabeculectomy. This is due to the fact that the level of blood flow before trabeculectomy was close to the baseline (healthy). In the hypertensive case the autoregulation is working properly only when IOP=28.7 mmHg, and the model predicts a value of blood flow that is higher in the HTwAR than in the HTwoAR. 51 On the other hand when IOP decreases, the model predicts an increase in blood flow over the baseline (healthy) value for both the HTwAR and HTwoAR cases. 52 3. COUPLED MODEL FOR VASCULAR REGULATION IN THE RETINA The retina exhibits vascular regulation (VR) of blood flow in response to changes in perfusion pressure and the metabolic needs of the tissue. Autoregulation (AR) is the ability to maintain relatively constant level of blood flow despite changes in perfusion pressure, whereas metabolic regulation (MR) is the ability to alter the level of blood flow in response to changes in tissue oxygen demand. AR and MR are achieved by varying the vascular tone of arterioles smooth muscles, causing the vessels to constrict or dilate. The model presented in Chapter 2 describes VR using a formula which phenomenologically regulates the level of blood flow in the retina. In this Chapter we improve this phenomenological approach by coupling the model described in Chapter 2 with a vessel wall mechanics model of blood flow VR in the retinal microcirculation [45]. This mechanistic model describes the response of vessels to local changes in pressure, shear stress, carbon dioxide and to the downstream metabolic state communicated via conducted responses. Although both models separately provide insight into key aspects of retinal blood flow regulation, the macrocirculation model is limited by using a phenomenological description of autoregulation, and the microcirculation model is limited by describing only the dynamics in the retinal microcirculation. In this Chapter, the microcirculation and macrocirculation models are combined into a single coupled model that includes a mechanistic description of blood flow VR in the microcirculation as well as the effects of IOP on the CRA and CRV. In Section 3.1 we describe the governing equations of the coupled model. In Section 3.2 we illustrate the solution algorithm for the coupled model. In Section 3.3 we present the model predictions of pressure, retinal blood flow and oxygen saturation 53 in response to changes in systemic blood pressure, intraocular pressure and tissue oxygen demand. In Section 3.4 we discuss the model results and limitations. The material and results presented in this Chapter are partially published in [118, 119] 3.1 Methods 3.1.1 Governing equations Retinal circulation is simulated using a mathematical model that describes the hemodynamic and mechanical properties of multiple vessel compartments, as depicted in Figure 3.1. As described in Chapter 2, the model is analogous to an electrical circuit in which blood is propelled through the system by a pressure gradient and the vessels are modeled as a network of resistors, representing the resistances to blood flow offered by the system. The model developed in this Chapter does not include capacitors since it focuses only on the steady state solution instead of on time-dependent dynamics. An analogy to Ohm’s Law is used to describe the blood flow (Q) through the system: Q = ∆P/R, where ∆P represents the pressure difference along a vessel and R is the vascular resistance. The model developed here is obtained by coupling two previously developed mathematical models for the retinal circulation. The first model, described in [45] (herein referred to as the microcirculation model), describes the blood flow through the large arterioles (LA), small arterioles (SA), capillaries (C), small venules (SV), and large venules (LV). Diameters of the LA and SA depend on the passive and active length-tension characteristic of the vessel wall. The active tension is generated by activation of the vascular smooth muscle according to changes in pressure, shear stress, oxygen (O2 ) metabolism and carbon dioxide (CO2 ). Thus, the microcirculation model provides a mechanistic description of retinal vascular regulation. The second model, described in Chapter 2 (herein referred to as the macrocirculation model), simulates blood flow upstream, within, and downstream of the microcir- 71 54 Microcirculation model SA LA P 2a P 1,2 R2a R1,IOP R1,LC P 2,3 R2b P 4a P 3,4 R3 LV SV R4a P 4,5 R4b intraocular region IOP P 1,post P 1,pre CRA R1 Rin C R5,IOP P 5,post LAMINA CRIBROSA RLTP P 5,pre CRV R5 P 1,in P 5,out Pin Pout P in R5,LC Rout P out translaminar region retrobulbar region Macrocirculation Model Figure 3.1. Schematic for retinal vascular network in the coupled model. The vasculature is divided into seven compartments: the central retinal artery (CRA), large arterioles (LA), small arterioles (SA), capillaries (C), small venules (SV), large venules (LV), and central retinal vein (CRV). Each compartment is modeled as a resistor (R). The retinal vasculature encompasses three regions: the intraocular region, translaminar region, and retrobulbar region. The vessels in the intraocular region are exposed Figure 3.1. toSchematic retinal vascular network the coupled model. intraocular for pressure (IOP), the vessels in theintranslaminar region are exposedistodivided an external that depends on the in the lamThe vasculature intopressure seven compartments: thestress central retinal ina cribrosa (shaded gray), the vessels in the retrobulbar region(C), are artery (CRA), large arterioles (LA),and small arterioles (SA), capillaries exposed to the retrolaminar tissue pressure (RLTp). Resistors marked small venules (SV), large venules (LV), and central retinal vein (CRV). with arrows indicate that their diameters change either passively (due to Each compartment is modeled a resistor (R). The retinal vasculature IOP) or actively (due toas vascular regulation). Dashed circled regions inencompasses three themodels intraocular region, translaminar region, dicate the regions: two previous (the macrocirculation model, presented in Chapter 2, and The the microcirculation [45]) that are combined in the and retrobulbar region. vessels in themodel intraocular region are exposed current work to form the coupled model. to intraocular pressure (IOP), the vessels in the translaminar region are exposed to an external pressure that depends on the stress in the lamina cribrosa (shaded gray), and the vessels in the retrobulbar region are exposed to the retrolaminar tissue pressure (RLTp). Resistors marked with arrows indicate that their diameters change either passively (due to IOP) or actively (due to vascular regulation). Dashed circled regions indicate the two previous models (the macrocirculation model, presented in Chapter 2, and the microcirculation model [45]) that are combined in the current work to form the coupled model. 55 culation. In particular, the macrocirculation model implements mechanical principles to describe blood flow through the CRA and the CRV. Coupling the models substitutes the phenomenological approach described in Section 2.1.6 with the mechanistic one used in [45], providing a more accurate description of retinal regulation in the context of a retinal vascular model that includes the hemodynamics in the vessels that supply (CRA) and drain (CRV) the retina. A complete description of the two original (i.e., uncoupled) models can be found in [45] (microcirculation model), and in Chapter 2 (macrocirculation model). In this Section, the relevant features of and changes to these two models necessary to generate the coupled model will be discussed. As in the previous two models, the coupled model is characterized by both variable and fixed resistances. The resistances Rin and Rout describe the blood vessel resistance upstream of the CRA and downstream of the CRV. Variable resistances are marked with an arrow in Figure 3.1, and the remaining resistances are fixed. Resistances R2a and R2b are vasoactive resistances that change according to VR mechanisms in the system. Resistances R1,LC , R1,IOP , R5,IOP and R5,LC are passive resistances that change according to the Law of Laplace (see Section 2.1.5) due to the external pressure of the IOP in the intraocular region and the mechanical stress exerted by the lamina cribrosa on the CRA and the CRV in the translaminar region (see Figure 3.1). The action of the lamina cribrosa on the translaminar portion of the CRA and CRV is described in Section 2.1.4 and in [48]. The coupled model differs from the macrocirculation model in the modeling assumption regarding the passive resistances: (i) the resistances R1,LC , R1,IOP , R5,IOP and R5,LC are modeled as compliant (not collapsible) tubes and the Starling effect is not included (see Section 2.1.5 for reference); (ii) venular resistances R4a and R4b are assumed to be constant. 56 Conservation of flow in the network yields the following nonlinear system of algebraic equations:     Pin − P1,post P1,post − P1,2   − =0   Rin + R1 + R1,LC (P ) R1,IOP (P )       P1,post − P1,2 P1,2 − P4,5   − =0  R1,IOP (P ) R2a (P , D, A) + R2b (P , D, A) + R3 + R4a + R4b   P1,2 − P4,5 P4,5 − P5,post    − =0   R2a (P , D, A) + R2b (P , D, A) + R3 + R4a + R4b R5,IOP (P )       P5,post − Pout   P4,5 − P5,post − = 0,  R5,IOP (P ) R1,LC (P ) + R5 + Rout (3.1) where the definitions of the pressures (P ) and resistances (R) are shown in Figure 3.1. The four variables of the system are given in the vector   P = P1,post , P1,2 , P4,5 , P5,post . (3.2) Pin and Pout are given data; the value of Pout is fixed at 14 mmHg, while the value of Pin is varied in the different simulations. Some of the vascular resistances depend on the pressure values (indicated in equation (3.1) by an explicit dependence in parentheses), while the remaining resistances are constant and equal to their control state values listed in Table 3.1. R2a and R2b also depend on the diameter (D) and smooth muscle activation (A) values obtained from system (3.3). 3.1.2 Modeling of vascular regulation Active changes in LA and SA resistances (R2a and R2b ) are achieved via four regulatory mechanisms that respond to changes in (a) systemic pressure (myogenic response), (b) wall shear stress (shear stress response), (c) adenosine triphosphate (ATP) concentration (metabolic response) and (d) carbon dioxide concentration (carbon dioxide response). 57 Table 3.1. Computed values of the vascular resistances at the control state. Value, Resistance Value, Resistance mmHg·s/mL mmHg·s/mL Rin 2.68 · 104 R4a 3.03 · 103 R1 8.60 · 103 R4b 1.86 · 103 R1,LC 1.97 · 102 R5,IOP 3.24 · 102 R1,IOP 1.01 · 103 R5,LC 6.09 · 101 R2a 8.21 · 103 R5 2.70 · 103 R2b 1.32 · 104 Rout 6.80 · 103 R3 6.63 · 103 The following equations are solved for the steady state diameter and smooth muscle activation of the LA and SA according to these four VR mechanisms:    dDi 1 Di   = Ti − Ttotal,i ,   ϕd T i  dt i = LA,SA, (3.3)     dAi 1   = Atotal,i − Ai ,  dt ϕa where Di and Ai are the diameter and level of smooth muscle activation, respectively. Constants ϕd = 1 s and ϕa = 1 min are the characteristic time constants regulating the changes in diameter and smooth muscle activation, and Di and T i are the control state values for diameter and wall tension. The wall tension Ti is given by the Law of Laplace as Ti = Pi − IOP Di , 2 i = LA,SA, (3.4) where Pi is the intraluminal pressure of the i-th compartment. The model representation for vessel wall tension Ttotal,i = Tpass,i + Atotal,i Tactive,i (3.5) 58 is the sum of a passive component (Tpass ) that arises from the structural elements of the vessel and a maximally active component (Tactive ) that arises from the constriction or relaxation of vascular smooth muscles. The passive component is assumed to be an exponential function of diameter [45] Tpass,i = Cpass,i eCpass,i [(Di /D0,i )−1] , 0 (3.6) and the maximally active component is modeled as a Gaussian function of vessel diameter [45] Tactive,i = Cact,i i2 h (D /D )−C 0 /C 00 e ( i 0,i act,i ) act,i . − (3.7) 0 0 00 The value of the constants Cpass , Cpass , Cact , Cact , Cact , and D0 , are taken directly from [45] and reported in Table 3.2. As defined in [45] the active component of tension in the vessel wall is the product of the maximally active and Atotal , which is a sigmoidal function (3.8) that varies between 0 and 1 and depends on a stimulus function (Stone , (3.9)) that determines the degree of vascular smooth muscles constriction: Atotal,i = 1 1 + e−Stone,i , i = LA,SA 00 Stone,i = Cmyo,i Ti − Cshear,i τi − Cmeta,i SCR,i − CCO2 ,i SCO2 ,i + Ctone,i (3.8) i = LA,SA. (3.9) The value of the constants Cmyo , Cshear , Cmeta , and CCO2 are obtained from [45] and 00 reported in Table 3.2. The constant Ctone is computed to guarantee that A = 0.5 in the control state. The stimulus function is a linear combination of the myogenic, shear-dependent, metabolic, and carbon dioxide VR mechanisms. • The myogenic mechanism contribution is Cmyo,i Ti where Ti represents the tension in the vessel wall computed as in equation (3.4). • The shear-dependent mechanism contribution is Cshear,i τi where τi is the wall shear stress τi = 32µi Qc,i , πDi3 (3.10) where µi represents the blood viscosity and Qc,i represents the blood flow in an individual vessel in compartment i. 59 Table 3.2. Values of the constant parameters for vascular regulation mechanisms. LA Parameter Cpass , dyn/cm SA Value Source 361.48 197.01 [45] 53.69 17.60 [45] 2114.2 3089.6 [45] 0 Cact 0.93 1.02 [45] 00 Cact 0.11 0.20 [45] Cmyo , cm/dyn 9.2 · 10−3 2.5 · 10−2 [45] Cshear , cm2 /dyn 2.58 · 10−2 2.58 · 10−2 [45] 200 200 [45] 8 · 10−4 1.31 · 10−4 [45] 0 Cpass Cact , dyn/cm Cmeta , (µM cm)−1 CCO2 , mmHg−1 00 Ctone D0 , µm 133.16*, 157.43** 60.18*, 70.54** 135.59 73.9 [45] * Values used for the results presented in Section 3.3.1 and in Section 3.3.2. ** Values used for the results presented in Section 3.3.3. 60 • The metabolic mechanism contribution is Cmeta,i SCR,i where SCR is a signal obtained by the integral of the ATP concentration cAT P (s) in plasma [45]. The ATP release depends on the oxygen saturation, SO2 of blood. tissue ζ tissue thickness, d vessel diameter, Di s vessel Figure 3.2. Krogh cylinder model schematic representation. Oxygen is assumed to diffuse from a vessel to a surrounding tissue cylinder of thickness d (left panel). A cross-sectional view of a group of adjacent parallel vessels is depicted on the right panel. • The carbon dioxide mechanism contribution is CCO2 ,i SCO2 ,i where SCO2 represents the carbon dioxide signal. This signal is a nonlinear function of the level of tissue CO2 content (cCO2 ) and the blood flow in the retinal circulation. More details on this mechanism can be found in [45]. Once the vessel diameter (D) is determined for the vessels in the LA and SA compartments, the values of the resistances are computed using Poiseuille’s Law R= 128µL , πnD4 (3.11) where L is vessel length, µ is blood viscosity, and n is the number of vessels in the compartment. In this model the transport of oxygen in the retinal vasculature and tissue is described using a Krogh cylinder type model. Oxygen is assumed to diffuse from a vessel to a surrounding cylinder of tissue, as in [45] (Figure 3.2). 61 3.1.3 Control state A control state represents conditions in a healthy eye. Control state values are indicated by an overbar. As in Section 2.1.2, the control state value of IOP is 15 mmHg, and the control state value of RLTp is 7 mmHg. The control state values of P in = 62.2 mmHg and P out = 14 mmHg are determined in the coupled model according to color Doppler imaging data (Section 2.1.2). To compute the control state value of Rin and Rout , the control state for the microcirculation model is computed first as in [45], using the value of P 1,2 and P 4,5 from Chapter 2, Table 2.4. Then, the pressure drop between the nodes P1,2 and P1,in and between the nodes P4,5 and P5,out is computed using Ohm’s Law. The value of Rin and Rout is equal to Rin = (P in − P 1,in )/Qun and Rout = (P 5,out − P 4,5 )/Qun , where all the quantities are equal to their control state value and Qun is the value of blood flow in the control state of the microcirculation model. The procedure used to compute the control state is iterative and is based on the control state assumptions described in Chapter 2 and [45]. First, the pressure values for P1,2 and P4,5 in Table 2.4 provide an initial guess for these pressure values, and then pressure drops across LA, SA, C, SV, and LV are computed as in [45]. The control state values of the LA and SA diameters are obtained by solving system (3.3), assuming that ALA = ASA = 0.5 in the control state. The control state value of the capillary diameter is assumed to be 6 µm, and a symmetry condition that dictates the length and number of vessels in each compartment is used to determine the control state diameter values of the SV and LV assuming that the wall shear stress values at control state are those reported in Table 3.3 [45]. The flow in each compartment is computed using Poiseuille’s Law, and then the numbers and lengths of the vessels in each compartment are computed [45]. A comprehensive list of the geometric and mechanical parameters characterizing the CRA and CRV are given in Chapter 2 and Table 3.3. 62 Table 3.3. Control state values for retinal vascular network compartments. Parameter CRA LA SA C SV LV CRV Retrobulbar 175 − − − − − 238 Translaminar 174.73 − − − − − 238.76 Retrobulbar 173.64 6 69.66 Diameter D, µm Length L, mm, total 105.54 48.07 155.70 235.18 10 5.86 4.27 0.54 4.27 5.86 10 Retrobulbar 8.8 − − − − − 8.8 Translaminar 0.2 − − − − − 0.2 Retrobulbar 10 5.86 4.27 0.54 4.27 5.86 10 1 4 38 191172 38 4 1 39.72 − − − − − 10.72 5.96 4.99 7.99 4.02 1.84 1.13 1.88 Retrobulbar 36.64 − − − − − 14.87 Translaminar 34.80 − − − − − 14.73 Retrobulbar 35.46 30 30 15 10 10 15.41 Retrobulbar 2.53 − − − − − 1.37 Translaminar 2.53 − − − − − 1.36 Retrobulbar 2.57 1.74 0.88 0.01 0.42 0.80 1.40 Number of segments n Wall thickness h, µm Pressure drop ∆P , mmHg Shear Stress τ , dyn/cm 2 Velocity v, cm/s Once the vascular network geometry for the model is determined, the control state resistance for each compartment is computed using Poiseuille’s Law (3.11). System (3.1) is then solved to compute the control state for the coupled model. The values of P1,2 and P4,5 for the coupled model are compared with the initial guess. The iterations cease when the difference between the initial guess and the computed value is within a certain tolerance. If the condition on the tolerance is not satisfied, the 63 predicted value is set as the new initial guess and the process is repeated. All control state values and parameters for the coupled model are listed in Tables 3.1, 3.3, 3.4. Table 3.4. Computed values of intraluminal blood pressures at the control state. Value, Pressure Value, Pressure mmHg 3.2 mmHg P in 62.22 P 3,4 22.97 P 1,in 45.94 P 4a 21.14 P 1,pre 40.71 P 4,5 20.01 P 1,post 40.60 P 5,post 19.81 P 1,2 39.98 P 5,pre 19.78 P 2a 34.99 P 5,out 18.13 P 2,3 27.00 P out 14.00 Solution algorithm In this Section we illustrate the solution procedure for coupling the macrocir- culation and microcirculation models (i.e. coupling the nonlinear problems in systems (3.1) and (3.3)). The problem is solved with the following fixed point algorithm: Let the values of IOP, RLTp, SP, DP be given, and let us set an initial guess for (i) the pressure at the nodes of the system in Figure 3.1, P 0 , (ii) the diameter and the 0 level of smooth muscle activation of the vessels in compartments LA and SA, DLA , 0 DSA , A0LA and A0SA , (iii) the value of oxygen saturation, ATP concentration and CO2 content at the inlet of the CRA, SO2 ,CRA (0), cAT P,CRA (0), cCO2 ,CRA (0), and (iv) the value of retinal blood flow Q0 . Let n = 0, . . . , N represent the iterations between the two systems and let us proceed as follows: 64 • Step 1: – Update Pin to be consistent with the given values of SP and DP – Solve thelamina cribrosa problem in [48] to obtain the value of PLC (ζ),  hLC hLC , with ζ ∈ − 2 2 For n ≥ 0 • Step 2: – Compute SOn 2 ,CRA (LCRA ), cnAT P,CRA (LCRA ), cnCO2 ,CRA (LCRA ) using [45] – Solve nonlinear system (3.3) for steady state solution (using Matlab ode n+1 n+1 n+1 solver) ⇒ DLA , DSA , An+1 LA and ASA n+1 n+1 – Compute the value of the resistances R2a , R2b using (3.11); • Step 3: n+1 n+1 n+1 n+1 – Solve nonlinear system (3.1) for P1,post P1,2 , P4,5 , and P5,post (using Mat- lab nonlinear solver) – Compute the value of retinal blood flow Qn+1 and the pressures at the nodes of the system in Figure 3.1 at iteration n + 1 using Ohm’s law; • Step 4: Check for convergence of the values of P1,2 and P4,5 ! |Pjn+1 − Pjn | < φ, max j={1,2},{4,5} |Pjn+1 | (3.12) where φ is a given tolerance; • Step 5: If (3.12) is satisfied then the quantities computed at step n + 1 are the solution of the coupled model. If (3.12) is not satisfied then set n = n + 1 and return to Step 2. 65 3.3 Results Here, we compare model predictions for two different versions of the model: the (uncoupled) microcirculation model and the coupled model defined in this Chapter. 1. In Sections 3.3.1 and 3.3.2 the coupled model and microcirculation model are used to predict the blood flow, blood oxygen saturation, and P1,2 and P4,5 values as the tissue oxygen demand M 0 = 1 cm3 O2 /100cm3 /min, and MAP and IOP vary in the range between (20, 165) mmHg and (7.5, 31) mmHg, respectively. A comparison between the results of the coupled model and microcirculation model is presented along with corresponding experimental measurements from rats [39, 120] and cats [121]. 2. In Section 3.3.3 the coupled model is used to predict blood oxygen saturation and blood flow as tissue oxygen demand varies in the range between (1, 5) cm3 O2 /100cm3 /min, and MAP= 93.3 mmHg, IOP= 15 mmHg and M 0 = 2.65 cm3 O2 /100cm3 /min. The model results are compared with clinical data on flicker stimulation and light-dark stimulation in which oxygen demand varies [30, 32, 34–36, 122]. 3.3.1 Effect of model coupling The effect of coupling the microcirculation and macrocirculation models is depicted in Figure 3.3 by comparing the microcirculation model (black curves) and coupled model (blue curves) predictions for P1,2 and P4,5 (panel a), the pressure drop along the CRA (panel b), and the total CRA resistance (panel c) as MAP is varied between 65 and 165 mmHg and IOP = 15 mmHg. A few key assumptions differ between the microcirculation model and the coupled model. For example, in the microcirculation model, P4,5 is held fixed (Figure 3.3a). The pressure drop along the CRA (i.e., ∆PCRA = Pin − P1,2 ) in the microcirculation model is not directly computed, but is assumed to be constant (Figure 3.3b), whereas the pressure drop along the CRA in 66 (a) A coupled model coupled model uncoupled model microcirculation model P1,2 60 50 40 30 P4,5 20 10 60 80 40 100 120 140 160 180 MAP (mmHg) coupled model 35 30 25 20 15 60 microcirculation uncoupled model model 80 100 120 140 160 180 MAP (mmHg) RCRAtot (mmHg s/mL) 70 C 45 Pin-P1,2 (mmHg) Pressure (mmHg) p 80 (b) B 90 x 10 5 (c) 4.5 4 microcirculation uncoupled model 3.5 model 3 2.5 2 1.5 coupled 1 model 0.5 0 60 80 100 120 140 160 180 MAP (mmHg) Figure 3.3. (a) Model predicted values of P1,2 (solid) and P4,5 (dashed) for the microcirculation model (black curves) and coupled model (blue curves). (b) Pressure drop along the CRA for the microcirculation model (black curve) and coupled model (blue curve). (c) Total resistance in the CRA in the microcirculation model (black curve) and coupled model (blue curve). In all panels, MAP is varied between 65 and 165 mmHg, IOP = 15 mmHg and tissue oxygen demand M0 = 1 cm3 O2 /100 cm3 /min. the coupled model varies according to the effects of the IOP-RLTp pressure difference and laminar effects on the CRA. Total resistance in the CRA is computed according to Ohm’s Law as RCRA = (Pin − P1,2 )/Q, where Q is the retinal blood flow through the system. Together, the panels in Figure 3.3 provide a mechanical rationale for why the model predictions for the microcirculation model and coupled model differ substantially, as shown in Figure 3.4. In Figure 3.4, the predictions for flow and venous oxygen saturation are compared for the two different models. In both cases, tissue oxygen demand is assumed to be held constant at 1 cm3 O2 /100 cm3 /min, MAP ranges from 65 to 165 mmHg, and IOP = 15 mmHg. In Figure 3.5, the various mechanisms of VR are turned “on” and “off” to demonstrate which mechanisms contribute most significantly to generating the VR plateau 67 A B (a) uncoupled microcirculation model model 3.5 venous saturation normalized retinal blood flow (b) 1 4 microcirculation uncoupled model model 3 2.5 2 1.5 coupled model 1 0.8 coupled model 0.6 0.4 0.2 0.5 0 60 80 100 120 140 MAP (mmHg) 160 180 0 60 80 100 120 140 MAP (mmHg) 160 180 Figure 3.4. (a) Predicted values of retinal blood flow (normalized) as MAP is varied (65 − 165 mmHg) using the microcirculation model (black curve) and coupled model (blue curve). (b) Predicted values for venous oxygen saturation as MAP is varied (65 − 165 mmHg) using the microcirculation model (black curve) and coupled model (blue curve). Arterial venous oxygen saturation was assumed to be 0.97 and tissue oxygen demand M0 = 1 cm3 O2 /100 cm3 /min. In both panels, IOP = 15 mmHg. in the coupled model. Consistent with the findings in [45], these results show that the metabolic and carbon dioxide mechanisms are most important for achieving VR. 3.3.2 Comparison with experimental measurements for changes in MAP and IOP To explore which version of the theoretical model (coupled or microcirculation) yields predictions that more closely mimic reality, predictions using both models are compared with experimental measurements from rats [39, 120] and cats [121]. In Figure 3.6, normalized values of blood flow are predicted in both the coupled model and microcirculation model, and compared with data from [39, 120, 121] as MAP is varied between 20 and 165 mmHg for IOP = 10 mmHg. 68 normalized retinal blood flow 2.5 myo+shear+meta+co2 myo+shear+meta myo+shear+co2 myo+shear 2 1.5 myo+shear +meta+co2 1 myo+shear 0.5 0 myo+shear+co2 myo+shear +meta 20 40 60 80 100 120 140 160 180 MAP (mmHg) Figure 3.5. Model predicted values of retinal blood flow (normalized) as MAP is varied (20 − 165 mmHg) using the coupled model. Vascular regulation mechanisms are turned on and off to assess their individual contributions. The following combinations are included: all mechanisms (blue), myogenic + shear + metabolic (green), myogenic + shear + CO2 (red), myogenic + shear (black curve). Here, IOP = 10 mmHg and tissue oxygen demand M0 = 1 cm3 O2 /100 cm3 /min. He et al. [39, 120] measure blood flow in rat retina at low-, normal- and elevatedblood pressure values corresponding to MAP = 59, 108, and 156 mmHg, respectively, as IOP is varied over a range of 10 − 110 mmHg. The coupled model and microcirculation model are only capable of handling IOP values between 7.5 and 31 mmHg, and these results are compared with the findings from [39, 120] in Figure 3.7. In the case of low MAP, as IOP increases, the diameter of the small arterioles decreases until the point of collapse, yielding the observed absence of blood flow when IOP > 22 mmHg (red, solid curve, Figure 3.7). 69 normalized retinal blood flow 4 data He data Tani coupled microcirc. uncoupled 3.5 3 microcirculation uncoupled model model 2.5 2 1.5 coupled model 1 0.5 0 20 40 60 80 100 120 140 160 180 MAP (mmHg) Figure 3.6. Normalized values of blood flow predicted by the coupled model (blue curve) and microcirculation model (black curve) are compared with data from [39,120] (red open markers) and [121] (magenta triangles) as MAP is varied between 20 and 165 mmHg for IOP = 10 mmHg and tissue oxygen demand M0 = 1 cm3 O2 /100 cm3 /min. 3.3.3 Comparison with experimental measurements for changes in oxygen demand Table 3.5 provides a summary of the observed effects of flicker stimulation and light-dark studies on retinal blood flow, blood oxygen saturation, and vessel diameter in multiple species (humans, rats, cats, and bullfrogs). These studies do not translate the degree of flicker stimulation or exposure to light into a quantifiable level of oxygen demand in the retina. In addition, the observed hemodynamic changes are not the same across species or for similar flicker stimulation frequencies or light exposure 70 Table 3.5. Summary of vascular response to flicker stimulation and light-dark adaptation studies. Study Species Results Flicker light stimula- humans - Retinal arterial and venous blood tion Source [30] flow increases during flicker stimulation - Major retinal arterial and venous diameter increases during flicker stimulation Flicker light stimula- humans tion - Retinal venous oxygen saturation [32] is higher during flicker stimulation than in light adaptation - Major retinal arterial and venous diameter increase during flicker stimulation Flicker light stimula- rats tion - Retinal venous oxygen saturation is [122] lower during flicker stimulation than in light adaptation Light-dark adaptation humans - Retinal arterial and venous oxygen [34] saturation is higher during dark than in light adaptation (in humans) - Oxygen consumption in the outer retina is higher in dark than light adaptation (in cats) - Oxygen consumption in the inner retina is similar in dark than light adaptation (in cats) Light-dark adaptation cats - Outer retina tissue oxygen partial [36] pressure is lower in dark than light adaptation (in cats) - Inner retina tissue oxygen partial pressure is higher in dark than light adaptation (in cats) - Retinal oxygen consumption is higher in dark than in light adaptation (in bullfrogs) Light-dark adaptation rats - Oxygen consumption in the outer and flicker light stimu- retina is higher in dark than in light lation adaptation - Inner retina activity is lower in dark than in light adaptation - Retinal blood flow is higher in light adaptation and during flicker stimulation than in dark adaptation [35] 71 4 coupled model 3.5 microcirculation uncoupled modelmodel normalized retinal blood flow data He 3 2.5 MAP = 156 mmHg 2 1.5 MAP = 108 mmHg 1 0.5 0 5 MAP = 59 mmHg 10 15 20 25 30 35 40 IOP (mmHg) Figure 3.7. Coupled model (solid) and microcirculation (dashed) model predictions of normalized blood flow as IOP is varied between 7.5 − 31 mmHg are compared with experimental measures [39,120] at low-, normaland elevated- blood pressure values corresponding to MAP = 59 (red), 108 (blue), and 156 mmHg (black), respectively. Note that these exact MAP values were provided by authors He et al. since average values were reported in [120]. All flow values are normalized with respect to the control state flow value for IOP = 15 mmHg, MAP = 93 mmHg and tissue oxygen demand M0 = 1 cm3 O2 /100 cm3 /min. time. Because of these aspects, these data provide mathematical modelers with a unique challenge of capturing the various responses according to the geometric and mechanistic assumptions built into their models. Figure 3.8 provides a first attempt of using a mathematical model to predict the change in perfusion when tissue oxygen demand is varied between 1 and 5 cm3 O2 /100 cm3 /min. The effects of different response mechanisms (myogenic, shear, metabolic, and carbon dioxide) are evaluated by running multiple model simulations with some 10 8 6 myo+shear myo+shear+meta myo+shear+CO2 Garhofer et al. 2004 (baseline) 4 1 1.5 2 2.5 3 3.5 4 4.5 5 3 3 oxygen consumption oxygen demand (cm O2/100 cm /min) 60 50 40 30 20 10 0 of the mechanisms active and some inactive. The model predicts asymmetric behavior in response to changes in tissue oxygen demand. For example, when all mechanisms are active, the model predicts a 6% decrease in perfusion when tissue oxygen demand is decreased by 50% and a 23% increase in perfusion when tissue oxygen demand is increased by 50%. If both the metabolic and carbon dioxide responses are impaired (light blue line), the model predicts a constant level of blood flow that does not respond to changes in oxygen demand, suggesting the importance of these two response mechanisms. ar h G Figure 3.8. Large arteriole (LA) blood flow vs. oxygen demand. The model predictions for different vascular regulation mechanisms are compared with clinical data [30] at baseline and during flicker stimulation. The following combinations are included: all mechanisms (blue curve), myogenic + shear + metabolic (black), myogenic + shear + CO2 (green), myogenic + shear (cyan). Here, IOP = 15 mmHg and oxygen consumption varies between 1 and 5 cm3 O2 /100 cm3 /min. C ha m ot e ta 2 70 l. myo+shear+ meta+CO2 12 Garhofer et al. 2004 (flicker) el 14 B m od l Large arteriole blood flow (mL/min) A % change in retinal blood flow during flicker 72 73 The predictions of blood flow through a large arteriole are also compared with data from a study on flicker stimulation in healthy individuals [30] in Figure 3.8. The clinical study measures: (i) the average value of blood flow in one large retinal arteriole of 7.4 ± 3.5 µl/min in light adapted conditions and (ii) a 59% increase in blood flow with flicker stimulation. In the model, a baseline value for oxygen demand is chosen to be M0 = 2.65 cm3 O2 /100 cm3 /min in order to yield the experimentally observed value of venous blood oxygen saturation of 60% [32]. The M0 value corresponding to flicker stimulation is obtained from a study showing that flicker stimulation induces a 48% increase in oxygen demand [122]. Using these clinical data points, Figure 3.8 shows that the model predictions are in good agreement with clinical data. The model predicted value of oxygen saturation in the large venules is also compared with data from a clinical study on flicker stimulation [32]. The clinical study measures an increase in venous saturation from 60%±5.7% to 64%±5.9% with flicker stimulation, while the model predictions show a decrease in venous saturation with flicker stimulation from 60% to 52%. 3.4 Discussion The vascular model of the retina developed in this study incorporates the me- chanical effects of IOP on the CRA and CRV (modeled as compressible tubes) and four VR response mechanisms to describe blood flow to the retina. The blood flow predicted by the coupled model is significantly closer to pressure-flow experimental data collected by [39, 120, 121] than blood flow levels predicted by the microcirculation model, supporting the improvements offered by the coupled model over the microcirculation model. This excellent agreement between the coupled model and experiments [39, 120, 121] provides confidence in the model assumptions and equations. 74 3.4.1 Effect of model coupling The coupled model accounts for the retinal macrocirculation and microcirculation. The coupling with the macrocirculation model affects the coupled model predictions which differ from the results obtained with the microcirculation model in the range of pressures where the VR plateau is absent. The coupled model and microcirculation model predict significantly different levels of intraluminal pressure directly upstream and downstream of the microcirculation compartments (i.e., P1,2 and P4,5 , as shown in Figure 3.3a). These differences result from differing model assumptions: in the microcirculation model P4,5 is held fixed, while in the coupled model, the passive effects of IOP and compression of the lamina cribrosa on the diameters of the CRA and CRV cause P1,2 and P4,5 to change. In addition, in the microcirculation model, the change in pressure along the CRA (i.e., ∆PCRA = Pin − P1,2 , as shown in Figure 3.3b) is not directly computed, but is held constant at 22 mmHg based on previous model predictions (see Chapter 2 and [45] for reference). To maintain this constant drop in pressure, the resistance in the CRA must decrease significantly (as shown in Figure 3.3c). In contrast, in the coupled model, the mechanical effects of IOP and the lamina cribrosa dictate ∆PCRA and cause nonlinear changes in the CRA resistance. As MAP increases, the change in pressure along the CRA also increases and the resistance of the CRA decreases slightly. The coupled model and microcirculation model were built to yield the same behavior in the control state, and thus the P1,2 and P4,5 values are the same for both models in the control state (as observed by the intersection of the curves at MAP = 93 mmHg in Figure 3.3a). As MAP increases above the control state, ∆PCRA is lower in the microcirculation model than the coupled model, and the reverse is true when MAP is decreased below the control state. Together, these relationships explain the differences in P1,2 and P4,5 predicted by the two models. Blood flow and oxygen saturation also differ between the coupled model and microcirculation model. While the two models predict nearly identical VR plateau (Fig- 75 ure 3.4a), the change in flow outside the pressure range for VR is more realistic in the coupled model, as evidenced by its agreement with experimental data (Figure 3.6). In the microcirculation model, flow is predicted to approach zero at a higher value of MAP than in the coupled model. Flow will equal zero once ∆P = P1,2 − P4,5 = 0. As can be seen in Figure3.3a, the curves for P1,2 and P4,5 in the microcirculation model intersect at a higher value of MAP than the coupled model, explaining the faster decline in flow as MAP decreases. The differences in flow outside the VR range for each model suggest that the assumptions in the microcirculation model are particularly relevant to study VR while the model assumptions in the coupled model provide a more accurate description of the overall hemodynamic behavior of the vascular supply to the retina. Similarly, the ranges of venous saturation predicted by both models are nearly equivalent (Figure 3.4b), but the coupled model yields a less steep decline in venous saturation as MAP decreases than the microcirculation model since blood flow does not decrease as rapidly. From the results in Figure 3.5, it is also apparent that the combined action of the metabolic and carbon dioxide responses is important in achieving flow regulation. 3.4.2 Comparison with experimental measurements for changes in MAP and IOP The comparison of the microcirculation and coupled model predictions with experimental data [39, 120, 121] highlights the importance of including the mechanical effects of IOP and the lamina cribrosa when modeling blood supply to the retina. As shown in Figure 3.6, the coupled model provides a very close match to experimental measures of flow in the retina of the rat [39, 120] and cat [121]. It is important to note that none of the data collected by [39, 120, 121] was used in the formulation or parameter optimization of this model; rather, the data set is only used to validate the model, and the agreement between theory and experiment highlights the value of including both a mechanical description of VR and IOP in the model. 76 Figure 3.7 provides further evidence that coupling the macrocirculation and microcirculation models is necessary to provide an accurate description of blood flow to the retina. Model simulations for IOP= 7.5−31 mmHg at low, normal, and elevated MAP values are compared with experimental data [39,120]. As both the coupled model and experimental data show, flow is relatively constant (with a slight decrease) over this range of IOP values. The microcirculation model predicts a much more significant decrease in flow with IOP which deviates substantially from the collected data. This decrease is due to the fixed levels of P1,2 and P4,5 (given a specific IOP value) in the microcirculation model, which leads to a decrease in flow as resistance increases with IOP. He et al. provide measures for IOP as high as 110 mmHg, but both theoretical models are limited once IOP is increased beyond a certain threshold (leading to the closure of the small arterioles in the coupled model). The coupled model is also unable to provide predictions for values of MAP or IOP that cause P4,5 = Pout . Of note, some work is still needed to fully account for the IOP effect (e.g. Starling resistor) on the retinal venules (SV and LV) and on the CRV. 3.4.3 Comparison with experimental measurements for changes in oxygen demand The increase in blood flow predicted by the model due to an increase in oxygen demand is not in the same proportion as the change in blood flow observed with the same decrease in oxygen demand, suggesting that vascular regulatory mechanisms may respond differently to different levels of oxygen demand. In addition, the model is used to quantify the relevance of the different regulatory mechanisms, suggesting that metabolic and carbon dioxide responses play a major role in achieving vascular regulation. While the model is able to quantify the average retinal response to changes in oxygen demand, it cannot differentiate between the contributions of the inner and the outer retina. Overall, the preliminary results of this model give important insight into metabolic flow regulation; however, the inability to capture all of the 77 in vivo observations listed in Table 3.5 [32, 34–36] suggests a need for more realistic descriptions of the retinal geometry and oxygen delivery to explain the response of blood flow to changes in retinal oxygen demand. To represent the retinal geometry more accurately, the model can be adapted to account for the multiple layers (vascular and avascular) of the retina. In this updated geometry, oxygen will be assumed to diffuse from retinal capillaries and the choroid into the various retinal tissue layers. The derivation of this extended model will be presented in Chapter 4. 78 4. MODELING OXYGENATION WITHIN MULTIPLE RETINAL LAYERS The model presented in this Chapter improves upon the models introduced in Chapters 2 and 3 by including more physiological representations of the following components: • retinal anatomy: the retinal tissue is defined using 6 different layers that include oxygen sources and sinks; • vascular architecture: capillaries are arranged in two layers at different depths providing oxygen to the tissue, together with the choroid (pure diffusion). • oxygen exchange: the exchange of oxygen occurs through the vessel walls due to the difference in oxygen concentration between blood and tissue. This boundary phenomenon requires a new formulation for the oxygen transport problem in the microcirculation (Section 4.2). In the next sections we will describe: (i) oxygen transport in the tissue (Section 4.1); (ii) oxygen transport in the microvasculature (Section 4.2); (iii) coupling between parts (i) and (ii) (Section 4.3); (iv) solution algorithm of the coupled problem (Section 4.4); and (v) preliminary results of numerical simulations (Section 4.5). 4.1 Oxygen transport in the retinal tissue The retinal tissue is described by the parallelepiped Ωt = Ωxy × (z0 , z6 ) in Fig- ure 4.1, where Ωxy represents the retinal surface, z0 = 0, and z6 = Lretina , where Lretina is the total thickness of the retina. The different retinal cells are grouped into 6 different layers according to the retinal anatomy, as depicted in Figure 4.1. Each 79 layer is labeled Ωt,j = Ωxy ×(zj−1 , zj ), with j = 1, . . . , 6. The values of the coordinates z0 , . . . , z6 and the measure of Ωxy are reported in Table 4.1. Layers from 1 − 3 comprise the inner retina and layers 4 − 6 comprise the outer retina. The thickness of the inner retina Lir and outer retina Lor are assumed to be Lretina equal, yielding Lir = Lor = [123]. 2 Layers 2, 4, 5, and 6 are avascular, whereas layers 1 and 3 contain the intermediate capillary plexus (referred to as C1 in the following sections) and the deep capillary plexus (referred to as C3 in the following sections), respectively (see Figure 4.1 for reference). x Ωxy VITREOUS y z0 LAYER 1 LAYER 2 inner nuclear layer and inner plexiform layer metabolically active z2 outer plexiform layer LAYER 3 deep capillary plexus (C3) metabolically active z3 Inner plexiform layer Inner nuclear layer Outer plexiform layer Outer nuclear layer inner segment of photoreceptors Membrana limitans externa Layer of rods and cones Pigmented layer z4 LAYER 5 Fibers of Müller outer nuclear layer LAYER 4 outer retina Membrana limitans interna Stratum opticum Ganglionic layer intermediate capillary plexus (C1) metabolically active z1 inner retina ganglion cell and nerve fiber layer metabolically active z5 LAYER 6 retinal pigmented epithelium and outer segment of photoreceptors z6 z CHOROID Figure 4.1. Left: Tissue layers of the retina clustered in 6 different layers in the Cartesian coordinate system (x, y, z). Right: Tissue layers of the retina, by Henry Vandyke Carter, Public Domain, via Wikimedia Commons, from [1], Plate 881. 80 The transport of oxygen in the retinal tissue is described by the stationary diffusion reaction equation in Cartesian coordinates (x, y, z) on Ωt , shown in Figure 4.1     −αt Dt ∆pt,j = Qj (pt,j , pc,j ), in Ωt,j , j = 1, . . . , 6     ∂pt,j   |z = 0, on Ωxy × {z = z0 }   ∂z 0     pt,j | − = pt,j | + , on Ωxy × {z = zj }, j = 1, . . . , 5 zj zj (4.1) ∂pt,j ∂pt,j   |−= | +, on Ωxy × {z = zj }, j = 1, . . . , 5   ∂z zj ∂z zj      pt,j |z6 = pt,choroid , on Ωxy × {z = z6 }       n · ∇ p = 0, on ∂Ω × (z , z ), xy t,j xy 0 6 where αt represents the oxygen solubility coefficient in the tissue, and pt,j represents the partial pressure of oxygen in the tissue region Ωt,j ; thus αt pt,j represents the ∂ is the partial oxygen concentration inside the tissue; ∆ is the Laplacian operator, ∂z derivative with respect to the z coordinate, n is the outward normal vector to ∂Ωxy × ∂ ∂ (z0 , z6 ) and ∇xy = ex + ey is the component of the gradient operator (∇) ∂x ∂y in the xy plane, where ex and ey are the unit vectors in the x and y direction, respectively. The quantity pc,j represents the partial pressure of free oxygen inside the capillaries in the capillary plexus Cj, Dt represents the oxygen diffusion coefficient in the tissue, and Qj (pt,j , pc,j ) represents the supply and consumption term of the equation. The boundary conditions in system (4.1) enforce: no oxygen flux at the vitreous level, continuity of partial pressure and flow at the internal interfaces between layers, Dirichlet condition at the level of the choroid, and no flux lateral condition on ∂Ωxy × (z0 , z6 ). The values of the constants in system (4.1) are reported in Table 4.1. The reaction term def Qj (pt,j , pc,j ) = Sj (pt,j , pc,j ) − Mj (pt,j ), (4.2) is the sum of (i) a source term Sj (pt,j , pc,j ) that represents the volumetric flow of oxygen supplied to the tissue by the capillaries in layers Cj, and (ii) a sink term Mj (pt,j ) that represents the consumption rate of oxygen due to the metabolic activity 81 of the tissue. Since the only layers containing capillaries are Layers 1 and 3, it follows that S2 = S4 = S5 = S6 = 0. (4.3) The explicit form of the source term Sj (j = 1, 3) will be discussed in Section 4.3. The consumption term Mj (pt,j ) is modeled by Michaelis–Menten kinetics [124] Mj (pt,j ) = αt Λmax pt,j j , pt,j + K1/2 (4.4) where Λmax represents the maximum rate of oxygen consumption of the tissue in j layer j, and K1/2 represents the partial pressure of oxygen at half of the maximum consumption. The maximum rate of oxygen consumption is assumed uniform in the inner retina, i.e. Λmax = Λmax = Λmax [124]. In the outer retina, the oxygen is 1 2 3 consumed primarily by the photoreceptors (Layer 5 in Figure 4.1). Consumption in = 0. The = Λmax Layers 4 and 6 is considered to be negligible [124, 125]; thus Λmax 6 4 control state values of the parameters Λmax and K1/2 are reported in Table 4.1. j Thus, the form of the reaction term in system (4.1) is:     S (p , p ) − Mj (pt,j ), j = 1, 3   j t,j c,j Qj (pt,j , pc,j ) = −Mj (pt,j ), j = 2, 5      0, j = 4, 6. 4.1.1 (4.5) Model reduction A model reduction of system (4.1) is preformed by averaging the problem on Ωxy , namely Z 1 αt Dt ∆pt,j dσ = Qj (pt,j , pc,j ) dσ, j = 1, . . . , 6. (4.6) |Ωxy | Ωxy Ωxy Z 1 def Here we will define pt,j (z) = pt,j dσ as the average value of pt,j at the retinal |Ωxy | Ωxy depth z, and split the Laplacian operator ∆ into its xy and z components, namely 1 − |Ωxy | Z ∆pt,j = ∆xy pt,j + ∂ 2 pt,j , ∂z 2 where ∆xy pt,j = ∂ 2 pt,j ∂ 2 pt,j + . ∂x2 ∂y 2 (4.7) 82 Table 4.1. Parameters for the oxygen transport in the retinal tissue Parameters Value Source Retinal surface |Ωxy |, cm2 10.9 [126] Retinal tissue thickness Lretina , µm 250 [124, 127] Inner retina thickness Lir , µm 125 [123] Outer retina thickness Lor , µm 125 [123] z1 − z0 , µm 25 [127] z2 − z1 , µm 75 [127] z3 − z2 , µm 25 [127] z4 − z3 , µm 67.5 [124] z5 − z4 , µm 25 [124] z6 − z5 , µm 32.5 [124] 2.4 · 10−5 [128] 2 · 10−5 [129] 80 [128] Layer 1 − 3, mmHg/s 26 [124] Layer 4, 6, mmHg/s 0 [124, 125] Layer 5, mmHg/s 90 [124] 2 [130] Thickness of retinal layers mLO2 Oxygen solubility coefficient in tissue αt , mL mmHg Diffusivity of oxygen in tissue Dt , cm2 /s Partial pressure of oxygen at the level of the choroid pt,choroid , mmHg Maximum rate of oxygen consumption Λmax j Partial pressure of oxygen at half maximum consumption K1/2 , mmHg The integral on left hand side (4.6) becomes (up to a negative sign): 1 |Ωxy | Z Z 1 ∂ 2 pt,j αt Dt ∆xy pt,j dσ + αt Dt dσ |Ωxy | Ωxy ∂z 2 Ωxy Z Z d2 1 1 αt Dt = ∆xy pt,j dσ + αt Dt 2 pt,j dσ |Ωxy | Ωxy dz |Ωxy | Ωxy Z d2 pt,j αt Dt ∇xy · (∇xy pt,j ) dσ + αt Dt = |Ωxy | Ωxy dz 2 Z  d2 pt,j 2 αt Dt    = n · ∇ p dσ + α D , j = 1, . . . , 6, xy t,j t t  |Ωxy | ∂Ω dz 2 xy  (4.8) 1 αt Dt ∆pt,j dσ = |Ωxy | Ωxy Z 83 where in step 1 the z component of the Laplacian operator can be moved outside of the integral since Ωxy is assumed uniform in z, and in step 2 the first term becomes an integral on ∂Ωxy using the divergence theorem. This term is zero due to the no lateral flux boundary condition in (4.1). The right hand side of equation (4.6) becomes: Z 1 Qj (pt,j , pc,j ) dσ = Qj (pt,j , pc,j ), |Ωxy | Ωxy j = 1, . . . , 6, (4.9) where Qj (pt,j , pc,j ) = S j (pt,j , pc,j ) − M j (pt,j ). The computation of the functional form of S j (j = 1, . . . , 6) will be discussed in Section 4.3 (equation (4.48)). The average oxygen consumption M is computed via a linear approximation as follows: Z pt,j Λmax Λmax pt,j (z) 1 j j αt dσ ' αt = Mj (pt,j ), M j (pt,j ) = |Ωxy | Ωxy pt,j + K1/2 pt,j (z) + K1/2 (4.10) for z ∈ (zj−1 , zj ), and j = 1, . . . , 6. The result of this averaging procedure is summarized in the following system which describes the transport of oxygen in the retinal tissue through the depth of the retina thickness (coordinate z):  d2 pt,j    −α D , = Qj (pt,j , pc,j ), in (zj−1 , zj ), j = 1, . . . , 6 t t   dz 2   dpt,j    | = 0, at z = z0   dz z0 pt,j |zj− = pt,j |zj+ , at z = zj , j = 1, . . . , 5     dpt,j dpt,j   |zj− = | +, at z = zj , j = 1, . . . , 5   dz dz zj     p |z = pt,choroid , at z = z6 . t,j 6 (4.11) The function pt,j = pt,j (z), thus the partial derivatives with respect to the z coordinate in (4.1) become ordinary derivatives. 4.2 Oxygen transport in the retinal microvasculature The retinal microcirculation model, presented in Chapter 3, is modified by split- ting the single capillary compartment in Figure 3.1 into two capillary compartments C1 and C3 located inside tissue Layers 1 and 3, respectively. It is assumed that 35% 84 (ψ1 = 0.35) of the total number of small arterioles and small venules branch into the intermediate capillary plexus (i.e., into tissue layer 1), and that 65% (ψ2 = 0.65) of the total number of small arterioles and small venules branch into the deep capillary plexus (i.e., into tissue layer 3). The ratio between the total number of capillaries in ψ1 compartments C1 and C3 is therefore equal to . To ensure similar hemodynamic ψ2 conditions in the two capillary compartments, the total blood flow is split using the same proportionality constant, namely: total blood flow comp. C1 = ψ1 ∗ (total retinal blood flow) (4.12) total blood flow comp. C3 = ψ2 ∗ (total retinal blood flow). A depiction of this updated microcirculation model is shown in Figure 4.2. The values of ψ1 and ψ2 are physiologically reasonable and similar to those used in [127], where ψ1 = 30%, ψ2 = 60% and the remaining 10% of arteries and veins were located at the level of the vitreous, and did not branch into capillaries. Here we do not account for the arterioles and veins that remain at the level of the vitreous. LA SA tissue Layer 1 VITREOUS SV LV C1 tissue Layer 2 CRA tissue Layer 3 C3 tissue Layers 4 − 6 Figure 4.2. Retinal microvasculature represented by the following compartments: the large arterioles (LA), small arterioles (SA), intermediate capillaries (C1), deep capillaries (C3), small venules (SV), and large venules (LV). The central retinal artery (CRA) is upstream of the model, and central retinal vein (CRV) is downstream of the model. CRV 85 Let us now focus on the transport of oxygen inside a vessel of the retinal microcirculation, such as that shown in Figure 4.3. The cylinder of length L and radius r with circular cross-section Σ in the polar coordinate system (ζ, θ, s) is used to represent a generic vessel in any of the microcirculation compartments depicted in Figure 4.2. r θ Σ s Σ ζ s=0 s=L Figure 4.3. Representative cylinder Ω in cylindrical coordinates (ζ, θ, s) with circular cross-section Σ, of length L and radius r. Oxygen transport is described by the stationary advection diffusion reaction equation in cylindrical coordinates on the domain Ω = [0, r) × [0, 2π) × (0, L) shown in Figure 4.3, namely ∇ · J = ϕ, in Ω, (4.13) where: def • J = uc − D∇c, is the oxygen flux; J = Jζ eζ + Jθ eθ + Js es , where eζ , eθ and es are the unit vectors in the ζ, θ and s (axial) direction, respectively; • u is the blood velocity, u = uζ eζ + uθ eθ + us es ; def • c = αc pc + HD c0 SO2 is the total oxygen concentration in the blood, which is the sum of (i) the oxygen dissolved in plasma and red blood cells (RBCs), denoted by αc pc , and (ii) the oxygen bound to hemoglobin in RBCs, denoted by HD c0 SO2 . The quantities pc and SO2 represent the partial pressure of free oxygen and the blood oxygen saturation, respectively. Oxygen saturation is related to pc through the Hill equation SO2 = pnH c , nH pc + pnH 50 (4.14) 86 where nH is the Hill coefficient and p50 is the partial pressure of free oxygen at 50% saturation. The constants αc , HD and c0 are the oxygen solubility coefficient, blood hematocrit and oxygen carrying capacity of the RBCs, respectively. Blood is modeled as a mixture of plasma and RBCs, therefore the solubility coefficient αc is computed as the weighted average of the plasma and RBC contributions: αc = HD α1 + (1 − HD )α2 , (4.15) where α1 represents the oxygen solubility in the RBCs and α2 represents the oxygen solubility in plasma. The values of αc , c0 , HD , α1 , α2 , p50 , and nH are reported in Table 4.2; • D is the diffusivity tensor; • ϕ is an oxygen sink or source in Ω. Table 4.2. Parameters for the oxygen transport in the retinal microvasculature Parameters Value Source Hematocrit HD 0.4 [131] Oxygen carrying capacity of RBCs c0 , mL O2 /mL mL O2 Oxygen solubility coefficient in RBCs α1 , mL mmHg mL O2 Oxygen solubility coefficient in plasma α2 , mL mmHg Diffusivity of oxygen in plasma Dζ , cm2 /s 0.5 [131] 3.38 · 10−5 [132] 2.82 · 10−5 [132] 2.18 · 10−5 [133] 1.5 · 10−7 [134] Partial pressure of free oxygen at 50% saturation p50 , mmHg 26 [131] Hill coefficient nH 2.7 [131] Diffusivity of hemoglobin in RBCs DHb , cm2 /s 87 Starting from (4.13), we will perform a model reduction using the following assumptions: • D is diagonal;   0  Dζ 0  D=  0 Dθ 0  0 0 Ds      (4.16) where Dζ , Dθ and Ds represent the diffusivity coefficients in the radial, angular and axial coordinate, respectively; • The right hand side of equation (4.13) is null, namely ϕ = 0, since no sources or sinks of oxygen are present in Ω. The exchange of oxygen, which filtrates through the vessel walls into the tissue, is modeled through the boundary conditions; • There is no dependence on the θ-coordinate ⇒ ∂· = 0, Jθ = 0, uθ = 0; ∂θ • A parabolic velocity profile à la Poiseuille is assumed ⇒ uζ = 0 ⇒ Jζ =  uζ c− ∂c , and us = 0 at ζ = r; Dζ ∂ζ • The contribution of hemoglobin radial diffusion in RBCs is neglected; !   DHb  ∂c ∂ Jζ = −Dζ = −Dζ α c pc + HD c0 SO2 , ∂ζ ∂ζ Dζ  (4.17)  where DHb represents the diffusivity of hemoglobin in RBCs. • The contribution of oxygen diffusion in the axial direction is negligible compared to the advective component ⇒   ((( D( ∂ Hb((( ( ( α( HD c0 SO2 , Js = us c − Ds (( c p( c+ ( Ds ((∂s (4.18) 88 The values of the constants Dζ and DHb are reported in Table 4.2. Given these assumptions, the term ∇ · J in equation (4.13) becomes:    1 ∂Jθ ∂Js 1 ∂ ζJζ + + ∇·J= ζ ∂ζ ζ ∂θ ∂s    1 ∂ ∂(αc pc ) ∂  =− ζDζ + us (αc pc + HD c0 SO2 ) . ζ ∂ζ ∂ζ ∂s (4.19) In this model, the exchange of oxygen between the blood vessel and the surrounding tissue is driven by the difference between the concentration of oxygen inside the vessel and the concentration of oxygen in the tissue surrounding the vessel. This physiological phenomenon is described by the wall-free model presented in [135] through the boundary condition of the system (4.20). Using the wall-free technique, the effect of the vessel wall is captured by the permeability coefficient Lp , yielding        1 ∂ ∂  ∂(α p ) c c   − ζD + u (α p + H c S ) = 0, in Ω  ζ s c c D 0 O 2   ζ ∂ζ ∂ζ ∂s    ∂(αc pc ) |ζ=r = Lp (αc pc |ζ=r − αt pt,w ), on Γ1 = {r} × [0, 2π) × (0, L) −Dζ   ∂ζ         on Σ0 = [0, r) × [0, 2π) × {0}.  pc = p̂c,0 , (4.20) The quantity pt,w represents the partial pressure of oxygen in the tissue surrounding the vessel, and p̂c,0 represents the partial pressure profile at the cylinder inlet. The parameter Lp regulates the exchange of oxygen between the blood inside the vessel and the surrounding tissue. The clinical literature does not provide a range of values for Lp ; thus in this work, the value of the parameter is estimated to yield physiologically reasonable oxygen concentration profiles in the mathematical simulations. The value of the parameter Lp is reported in Table 4.3. 4.2.1 Model reduction A model reduction of the system (4.20) is performed by integrating over the crosssection Σ = [0, r) × [0, 2π). Further assumptions on us , pc and SO2 are necessary. It 89 is assumed that the dependence on ζ and s for each of the quantities of interest can be expressed as the product of (i) the quantity cross-sectional average (overlined function) that depends on the axial coordinate s only, and (ii) a shape function f that depends on the radial component ζ only. So we write: Z 2π Z r 1 us ζ dζ dθ us (ζ, s) = us (s)fu (ζ), where us (s) = 2 πr 0 0 Z 2π Z r 1 pc (ζ, s) = pc (s)fp (ζ), where pc (s) = 2 pc ζ dζ dθ πr 0 0 Z 2π Z r 1 SO2 (ζ, s) = S O2 (s)fSO2 (ζ), where S O2 (s) = 2 SO2 ζ dζ dθ. πr 0 0 (4.21) The shape function can be computed explicitly by enforcing compatibility with the boundary conditions. The assumption of Poiseuille flow for the blood implies i Parabolic profile of us fu (ζ) = aζ 2 + bζ + c. ⇒ ⇒ ii The maximum value of us is attained at ζ = 0 fu0 (ζ)|0 = (2aζ + b)|0 = 0 ⇒ iii fu (ζ)|r = 0 ar2 + c = 0 ⇒ ⇒ fu0 (ζ)|0 = 0 b = 0. (4.22) c = −ar2 . The separation of variables in equation (4.21) implies that fu (ζ) must satisfy Z 2π Z r fu (ζ)ζ dζ dθ = πr2 , (4.23) 0 0 which implies that, πr2 = Z 2π Z r Z 2π Z fu (ζ)ζ dζ dθ = 0 0 0 2 Therefore fu (ζ) = 2 1 − 2 a = − 2, r ! 0 r " ζ 4 ζ 2 r2 a(ζ 2 − r2 )ζ dζ dθ = 2πa − 4 2 ⇒ #r 0 c = 2. (4.24) ζ , which corresponds to the classic Poiseuille profile. r2 90 The boundary condition for the partial pressure of oxygen pc in system (4.20) requires additional attention. Replacing pc in the second condition of system (4.20) with the expression in (4.21) gives: −Dc ∂(αc pc ) |Γ1 = Lp (αc pc |Γ1 − αt pt,w ) ∂ζ −Dc ∂(αc pc (s)fp (ζ)) |Γ1 = Lp (αc pc (s)fp (ζ)|Γ1 − αt pt,w ) ∂ζ (4.25) −Dc αc pc (s)fp0 (r) = Lp (αc pc (s)fp (r) − αt pt,w ). This is a differential equation for fp (r) that can be solved if the dependence on the coordinate s is eliminated. To this end we define the function p∗ (ζ, s) = pc (ζ, s) − βpt,w (s), def where β = αt , αc (4.26) representing an extension of the tissue oxygen partial pressure pt,w from the vessel wall into Ω [127, 136]. Note that pt,w is assumed to be constant on each cross-section and so it does not depend on ζ. Consistent with the equations (4.21), it is assumed that p∗ (ζ, s) = p∗ (s)f ∗ (ζ), and so pc (s)fp (ζ) = p∗ (s)f ∗ (ζ)+βpt,w (s). Moreover, since (4.27) Z 0 2π Z r f ∗ (ζ)ζ dζ dθ = πr2 0 we can write 1 pc (s) = 2 πr = Z 0 2π Z 0 r 1 pc (ζ, s)ζ dζ dθ = 2 πr Z 0 2π Z r  ∗  p (s)f ∗ (ζ) + βpt,w (s) ζ dζ dθ 0  1 ∗ 2 2 p (s)πr + βp (s)πr = p∗ (s) + βpt,w (s). t,w πr2 k ! (4.28) k+2 ζ 1 − δ k , where k + 2(1 − δ) r ∗ k and δ are arbitrary constants [136] that can be chosen so that f (ζ) satisfies the Consider the candidate shape function f ∗ (ζ) = second equation in (4.20). We verify the compatibility of f ∗ (ζ) with the second equation in (4.20) and with condition (4.23) in the following way: 91 i Z 0 2π Z 0 r ! k k + 2 ζ f ∗ (ζ)ζ dζ dθ = 1 − δ k ζ dζ dθ r 0 k + 2(1 − δ) 0 " #r ζ2 δ ζ k+2 k+2 − k = 2π k + 2(1 − δ) 2 r k+2 !0 k+2 2 k+2 r δ r = 2π − k k + 2(1 − δ) 2 r k+2 Z 2π = 2 π = πr Z r (4.29) 2r2 + kr2 − 2δr2   k + 2(1 − δ) 2(k + 2)   k + 2 ( ((( −( δ) (. (( k( +( 2(1 −( δ) ( k( + 2(1 2( ii Assuming a parabolic profile for f ∗ (ζ), we set k = 2. Thus ! 2 2 ζ 2(1 − δ) f ∗ (ζ) = 1−δ 2 , f ∗ (r) = , 2−δ r 2−δ  ∗ 0  ∗ 0 4δ 4δζ , f (r) = − , f (ζ) = − (2 − δ)r2 (2 − δ)r (4.30) iii The last step is to determine the value of δ such that the second equation in (4.20) is verified: −Dζ ∂(αc pc ) |Γ1 = Lp (αc pc |Γ1 − αt pt.w ) ∂ζ −Dζ ∂(αc pc ) |Γ1 = Lp (αc pc − αt pt.w )|Γ1 ∂ζ − αc Dζ ∂(p∗ + βpt.w ) |Γ1 =  αc Lp (p∗ )|Γ1 ∂ζ pt,w =pt,w (s) *  ∂p∗ ∂βp t.w −Dζ |Γ1 + D |Γ1 = Lp (p∗ )|Γ1  ζ ∂ζ ∂ζ  −Dζ  ∗ (ζ)) ∂( p∗ (s)f  ∗ (ζ))| |Γ1 = Lp ( p∗ (s)f Γ1 ∂ζ −Dζ [f ∗ ]0 (r) = Lp f ∗ (r) δ= Lp r 2Dζ + Lp r ⇒ (δ 6= 2). 4δDζ 2(1 − δ) = Lp (2 − δ)r 2−δ (4.31) 92 ! 2 2 ζ The function f ∗ (ζ) = 1 − δ 2 , with δ defined by (4.31) is compatible 2−δ r with the wall-free boundary condition and satisfies equation (4.23). Therefore this function is chosen to be the shape function for p∗ (ζ, s) in (4.27). Oxygen saturation (SO2 ) is computed using the Hill equation as a function of the free oxygen partial pressure (pc ) and it is related to the oxygen bound to hemoglobin in the RBCs. This current work does not focus on the spatial distribution of RBCs inside of the vessel lumen; thus oxygen saturation is assumed to be a function of pc and fSO2 (ζ) ≡ 1. In this way, condition (4.23) is automatically verified by fSO2 (ζ), and def SO2 (ζ, s) = S O2 (s) = pnH c (s) . nH pc (s) + pnH 50 (4.32) Now let’s substitute p∗ into system (4.20) and then proceed to the model reduction       1 ∂ i ∂(αc (p∗ + βpt,w )) ∂ h  ∗  − u ζD + α (p + βp ) + H c S = 0, in Ω  s ζ c t,w D 0 O 2   ζ ∂ζ ∂ζ ∂s    ∂(αc (p∗ + βpt,w )) = Lp αc p∗ , on Γ1 −D ζ   ∂ζ        ∗  on Σ0 .  p = p̂c,0 − βpt,w , (4.33) The cross-sectional integration of the radial diffusive contribution in (4.33) gives    r Z 2π Z r 1 ∂ ∂(αc (p∗ + βpt,w )) ∂(αc (p∗ + βpt,w )) − ζDζ ζ dζ dθ = −2π ζDζ ∂ζ ∂ζ 0 0 ζ ∂ζ 0 ∂(αc (p∗ + βpt,w )) |Γ1 = −2πrDζ ∂ζ = 2πrαc Lp p∗ |Γ1 = 2πrαc Lp p∗ (s)f ∗ (r), (4.34) and the cross-sectional integration of the axial contribution in (4.33) gives Z 2π Z r i ∂ h us αc (p∗ + βpt,w ) + HD c0 SO2 ζ dζ dθ. 0 0 ∂s The integral in (4.35) is the sum of three contributions: (4.35) 93 1. 2π Z r Z 0 0 ∂ d [αc us p∗ ] ζ dζ dθ = ∂s ds Z 2π Z 0 r   us (s)fu (ζ)αc p∗ (s)f ∗ (ζ) ζ dζ dθ 0  d = αc us (s)p∗ (s) ds = αc Z 2π Z 0 r   fu (ζ)f ∗ (ζ) ζ dζ dθ 0  d us (s)p∗ (s) πr2 ωu∗ , ds (4.36) where 1 ωu∗ = 2 πr Z 2π Z r  fu (ζ)f ∗ (ζ) ζ dζ dθ 0 0  ! ! Z 2π Z r 2 2 2 ζ 1 2 1 − ζ 1 − δ 2  ζ dζ dθ = 2 2 πr 0 r 2−δ r 0 ! Z r 8 π ζ5 ζ3 = δ 4 − (1 + δ) 2 + ζ dζ (2 − δ) πr2 0 r r " #r 8 ζ6 ζ4 ζ2 = δ − (1 + δ) 2 + (2 − δ)r2 6r4 4r 2 0   8 2δ − 3 − 3δ + 6 6 2(3 − δ) = . r = (2 − δ)r2 12r4 3(2 − δ) def  (4.37) 2. Z 0 2π Z 0 r  ∂  d αc βus pt,w ζ dζ dθ = ∂s ds Z 0 2π Z r   αc βus (s)fu (ζ)pt,w (s) ζ dζ dθ 0  d = αc β us (s)pt,w (s) ds = αc β Z 2π 0 Z r   fu (ζ) ζ dζ dθ 0  d us (s)pt,w (s) πr2 . ds (4.38) 94 3. Z 0 2π Z 0 r ∂ d [us HD c0 SO2 ] ζ dζ dθ = ∂s ds Z 0 2π Z rh i HD c0 us (s)fu (ζ)S O2 (s) ζ dζ dθ 0  Z 2π Z r   d  = HD c0 us (s)S O2 (s) fu (ζ) ζ dζ dθ ds 0 0 = HD c0  d us (s)SO2 (s) πr2 . ds (4.39) Notice that in (4.36), (4.38), and (4.39) the partial derivative with respect to s can be moved outside of the double integral since the cross-section Σ does not depend on the axial coordinate s. In addition, outside of the integral the partial derivative becomes a total derivative with respect to s, since the double integral removes the dependence from the radial coordinate ζ. The sum of (4.34), (4.36), (4.38), and (4.39) gives:    d 2 ∗ us (s)πr αc p (s)ωu∗ + αc βpt,w (s) + HD c0 S O2 (s) + 2πrαc Lp p∗ (s)f ∗ (r) = 0, ds (4.40) where the second equation in (4.20) is now included in (4.40) as a reaction term. The term us (s)πr2 = Qc (s) is equal to the blood flow in the vessel Ω. Note that Qc is determined by the vascular resistance and the pressure drop along the vessel (see Ohm’s law in Chapter 2 for reference), and so it does not depend on the axial coordinate s. Thus, equation (4.40) is written as:   d  ∗ Qc αc p (s)ωu∗ + αc βpt,w (s) + HD c0 S O2 (s) + 2πrαc Lp p∗ (s)f ∗ (r) = 0. (4.41) ds To complete the model reduction procedure, equation (4.40) is written in terms of the original variable pc :   d    Q α p (s)ω + H c S (s) +  c c c u∗ D 0 O2  ds     d  +Qc αc β (1 − ωu∗ ) pt,w + ds     +2πrαc Lp pc (s) − βpt,w f ∗ (r) = 0,      p (0) = p , c c,0 (4.42) for s ∈ (0, L) 95 where pc,0 is the value of pc when s = 0. System (4.42) describes the reduced one dimensional problem for the transport of oxygen in one vessel of the retinal microcirculation depicted in Figure 4.2 and 4.3. 4.3 Coupling between retinal tissue and retinal microvasculature To describe the oxygenation of the retinal microvasculature and the retinal tissue, the systems are coupled using the following procedure: The first equation in system (4.42) is written for every compartment in Figure 4.2: Qc,j  d  αc pc,j (s)ωu∗,j + HD c0 S O2 ,j (s) + ds   d  +Qc,j αc β 1 − ωu∗,j pt,w,j + ds   +2πrj αc Lp,j pc,j (s) − βpt,w,j fj∗ (rj ) = 0, for s ∈ (0, Lj ), (4.43) j = LA,SA,C1,C3,SV,LV, where the subscript j denotes the compartment where each quantity is computed. The values of the quantities introduced in equation (4.43) are reported in Table 4.3. Table 4.3. Geometrical and physiological parameter of the network in Figure 4.2 Parameters VALUE Source LA SA C1 C3 SV LV radius r, µm 52.8 24 3 3 34.8 77.8 ch. 3 Length L, cm 0.59 0.43 0.054 0.054 0.43 0.59 ch. 3 4 38 66910 124262 38 4 ch. 3 0.5 1.1 21.6 21.6 0 0 - Number of vessels N Permeability coefficient Lp , 10−4 cm/s 96 The quantity pt,w,j in equation (4.43) represents the partial pressure of oxygen in the tissue surrounding the vessels in compartment j and is computed as follows: pt,w,j = pt,1 (z = 0) for j = LA,SA,SV,LV ∼ pt,w,j = pt,j for j = C1,C3 (4.44) where: ∼ pt,C1 ∼ pt,C3 1 def = z1 − z0 Z 1 = z3 − z2 Z def z1 pt,1 (z) dz z0 z3 pt,3 (z) dz. z2 For compartments LA, SA, SV, and LV, the quantity pt,w,j is set equal to the value of partial pressure of oxygen at the level of the vitreous, while in compartments C1 and C3, the quantity pt,w,j is equal to the average of pt,j computed along the tissue z coordinate, as shown in (4.44). The modeling choices for pt,w,j , summarized in (4.44), imply that the quantity does not depend on the axial coordinate s, so dpt,w,j = 0, ds j = LA,SA,C1,C3,SV,LV, (4.45) The last term of equation (4.43) represents an oxygen volumetric rate per unit of length and is defines as:   qj (s) = 2πrj αc Lp,j pc,j (s) − βpt,w,j fj∗ (rj ), def for j = LA,SA,C1,C3,SV,LV. (4.46) Following the same assumptions of the model described in Chapter 3, it is assumed that compartments SV and LV do not exchange oxygen with the surrounding tissue. This assumption is implemented by setting Lp,SV = 0 and Lp,LV = 0, as reported in Table 4.3. This implies that qSV (s) ≡ 0 and qLV (s) ≡ 0. 97 Combining (4.43), (4.44), (4.45), and (4.46) into system (4.42) gives:   d     Qc,j αc pc,j (s)ωu∗,j + HD c0 S O2 ,j (s) = −qj (s), s ∈ (0, Lj ],   ds     pc,LA (0) = pc,0 ,         NLA Qc,LA cLA (LLA ) = NSA Qc,SA cSA (0), NSA Qc,SA cSA (LSA ) = NC1 Qc,C1 cC1 (0) + NC3 Qc,C3 cC3 (0),      cC1 (0) = cC3 (0),       NC1 Qc,C1 cC1 (LC1 ) + NC3 Qc,C3 cC3 (LC1 ) = NSV Qc,SV cSV (0),       NLA Qc,SV cSV (LSV ) = NLV Qc,LV cLV (0), (4.47) where cj (s) = αc pc,j (s)ωu∗,j + HD c0 S O2 ,j (s), and Nj is the number of parallel vessels in the compartment j, with j = LA,SA,C1,C3,SV,LV. The last five equations in system (4.47) guarantee the conservation of oxygen between the different compartments. To conclude the coupling procedure between the retinal tissue and the retinal microcirculation, we provide the explicit formulation of the sources of oxygen in the j-th retinal layer, S j in (4.11) where j = 1, 3, Z LC1 Z LC3 NC1 NC3 def def S1 = qC1 (s) ds, S3 = qC3 (s) ds. (z1 − z0 ) |Ωxy | 0 (z3 − z2 ) |Ωxy | 0 (4.48) 4.4 Solution algorithm Let us define the linearized quantities M̂ and ŜO2 . The nonlinear consumption term M i (pt,i ) (i = 1, . . . , 6) in equation (4.10) is lin- earized as follows: Mi (pt,i ) = αt Λmax pt,i (z) αt Λmax i i = p (z) = M̂i (pt,i ) pt,i (z), pt,i (z) + K1/2 pt,i (z) + K1/2 t,i (4.49) αt Λmax i , and z ∈ (zj−1 , zj ), j = 1, . . . , 6. pt,i (z) + K1/2 The oxygen saturation S O2 ,j (s) (j = LA,SA,C1,C3,SV,LV) in (4.42) is linearized def where M̂i (pt,i ) = as follows: pnH−1 pnH c,j c,j S O2 ,j (s) = nH = nH pc,j = ŜO2 ,j (pc,j )pc,j (s), nH pc,j + p50 pc,j + pnH 50 (4.50) 98 pnH−1 c,j . nH pnH c,j + p50 Consider the coupled problem consisting of the nonlinear systems (4.11) and (4.47). def where ŜO2 ,j (pc,j ) = The coupled problem is solved following the algorithm depicted in Figure 4.4. The algorithm consists of two nested fixed point iterations. The outer loop (index k) enforces the conservation of mass between the two systems while the inner loops (indexes n and m) solve the nonlinear problems of the transport of oxygen in the microcirculation and in the tissue systems. The algorithm proceed as follows: Given the values of Qc,j and ωu∗,j , and the initial guess for p0c,j (s) s ∈ [0, Lj ] (j = LA, . . . , LV) and p0t,i (z) z ∈ [zi−1 , zi ] (i = 1, . . . , 6): 0 • Step 1: Compute qj0 (s) (j = LA, . . . , LV) and S i (i = 1, 3) with formula (4.46) and (4.48), respectively; for k ≥ 0 • Step 2: k – Set n = 1 and pk,0 c,j (s) = pc,j (s) , (j = LA, . . . , LV) k – Set m = 1 and pk,0 t,i (s) = pt,i (s) , (i = 1, . . . , 6); for n ≥ 1 and m ≥ 1 • Step 3: (j = LA, . . . , LV) k,n−1 − Set pk,n (s) c,j (s) = pc,j  nH−1 pk,n (s) c,j (s) =  ; − Compute ŜOk,n nH 2 ,j k,n nH pc,j (s) + p50 99 k=0   s ∈ 0, Lj initialize: given: p0c,j (s) p0t,i (z) z ∈ [zi−1 , zi ]   compute: qjo (s) s ∈ 0, Lj 0 Si 0 n = 1, pk,0 c,j = pc,j 0 m = 1, pk,0 t,i = pt,i oxygen transport oxygen transport in microcirculation in retinal tissue k,n−1 set pk,n c,j = pc,j k,m−1 set pk,m t,i = pt,i compute ŜOk,n 2j compute M̂ik,n solve (4.51) ⇒ pk,n+1 c,j solve (4.53) ⇒ pk,m+1 t,i n=n+1 m=m+1 e(pk,n+1 , pk,n c,j c,j ) < φc,j no e(pk,m+1 , pk,m t,i t,i ) < φt,i no update: yes compute pk+1 c,j yes compute pk+1 t,i compute qjk+1 k+1 compute S i k =k+1 k =k+1 n=1 set pk,0 c,j = m=1 e(qjk+1 , qjk ) < φq,j pkc,j k+1 e(S i no k set pk,0 t,i = pt,i k , S i ) < φS,i k e(pk+1 c,j , pc,j ) < φc,j k e(pk+1 t,i , pt,i ) no < φt,i yes stop: pc,j = pk+1 c,j pt,i = pk+1 t,i Figure 4.4. Schematic of the solution algorithm for the coupled problem of oxygen transport. In the diagram j =LA,SA,C1,C3,SV,LV and i = 1, . . . , 6. 100 • Step 4: Solve the linearized version of system (4.47) for pk,n+1 (s) c,j   d  k,n+1  k,n k,n+1   Q α p (s)ω + H c Ŝ (s)p (s) = −qjk (s), s ∈ (0, Lj ], c,j c c,j u∗,j D 0 O2 ,j c,j   ds   k,n   pk,n+1  c,LA (0) = pc,LA (0),    k,n+1 k,n+1     NLA Qc,LA cLA (LLA ) = NSA Qc,SA cSA (0), k,n+1 NSA Qc,SA ck,n+1 (0) + NC3 Qc,C3 ck,n+1 (0), SA (LSA ) = NC1 Qc,C1 cC1 C3      ck,n+1 (0) = cC3 (0)k,n+1 ,  C1      NC1 Qc,C1 ck,n+1 (LC1 ) + NC3 Qc,C3 ck,n+1 (LC1 ) = NSV Qc,SV ck,n+1 (0),  C1 C3 SV      NLA Qc,SV ck,n+1 (LSV ) = NLV Qc,LV ck,n+1 (0), SV LV (4.51) (s) = αc pk,n+1 (s)ωu∗,j +HD c0 ŜOk,n (s)pk,n+1 (s), and j = LA, . . . , LV; where ck,n+1 j c,j c,j 2 ,j • Step 5: Check convergence of the solution of (4.51): " # k,n+1 k,n max |p (s) − (s)| p s∈[0,L ] c,j c,j j e(pk,n+1 , pk,n max < φc , (4.52) c c ) = k,n+1 j=LA,...,LV maxs∈[0,Lj ] |pc,j (s)| where φ is a given tolerance; • Step 6: If (4.52) is satisfied go to Step 7. If (4.52) is not satisfied, then set n = n + 1 and return to Step 3; • Step 7: (i = 1, . . . , 6) k,n−1 − Set pk,n (z) t,i (z) = pt,i − Compute M̂jk,n (z) = αt Λmax i pk,n t,i (z) + K1/2 ; • Step 8: Solve the linearized version of system (4.11) for pk,n+1 (z) t,i    d2 pk,n+1 k t,i   −αt Dt , = S i − M̂ik,n pk,n+1 , in (zj−1 , zj ), i = 1, . . . , 6  t,i 2  dz     dpk,n+1 t,i   at z = z0   dz |z0 = 0, pk,n+1 |zi− = pk,n+1 |zi+ , at z = zi , i = 1, . . . , 5  t,i t,i   k,n+1 k,n+1   dpt,i dpt,i   |zi− = |zi+ , at z = zi , i = 1, . . . , 5   dz dz      pk,n+1 |z6 = pt,choroid , at z = z6 ; t,i (4.53) 101 • Step 9: Check convergence of the solution of (4.53): " # k,n+1 k,n |p max (z) − p (z)| z∈[zi−1 ,zi ] t,i t,i , pk,n < φt e(pk,n+1 t t ) = max k,n+1 i=1,...,6 maxz∈[zi−1 ,zi ] |pt,i (z)| (4.54) where φ is a given tolerance; • Step 10: If (4.54) is satisfied go to Step 11. If (4.54) is not satisfied, then set m = m + 1 and return to Step 7; • Step 11: – Update the variables in the fixed point iterations on the index k using the relaxation method with parameters λc and λt ∈ (0, 1]: k,n+1 pk+1 (s) + (1 − λc ) pkc,j (s), c,j (s) = λc pc,j s ∈ [0, Lj ], k,n+1 pk+1 (z) + (1 − λt ) pkt,i (z), t,i (z) = λt pt,i z ∈ [zi,1 , zi ], k+1/2 – Compute qj k+1/2 (s) and S i j = LA, . . . , LV, i = 1, . . . , 6 (4.55) with equations (4.46) and (4.48), respec- k+1 tively, using the values pk+1 c,j and pt,i – Update q and S using the relaxation method with parameters λq and λS ∈ (0, 1] k+1/2 qjk+1 (s) = λq qj k+1 Si = k+1/2 λS S i  (s) + 1 − λq qjk (s), k λS ) S i , + (1 − s ∈ [0, Lj ], j = LA, . . . , LV, i = 1, 3; (4.56) • Step 12: Check convergence: " e(pk+1 , pkc ) = c max k maxs∈[0,Lj ] |pk+1 c,j (s) − pc,j (s)| maxs∈[0,Lj ] |pk+1 c,j (s)| j=LA,...,LV " e(pk+1 , pkt ) = max t k maxz∈[zi−1 ,zi ] |pk+1 t,i (z) − pt,i (z)| " max j=LA,...,LV maxs∈[0,Lj ] |qjk+1 (s) − qjk (s)| maxs∈[0,Lj ] |qjk+1 (s)| < φc (4.57) # maxz∈[zi−1 ,zi ] |pk+1 t,i (z)| i=1,...,6 e(q k+1 , q k ) = # < φt (4.58) < φq (4.59) # 102  e(S k+1 k ,S ) =  k+1 k |S −S | max  i k+1 i  i=1,3 |S i | < φS ; (4.60) • Step 13: If (4.57), (4.58), (4.59), and (4.60) are verified then pk+1 c,j (s) and pk+1 t,i (z), for s ∈ [0, Lj ], j = LA, . . . , LV, z ∈ [zi−1 , zi ], i = 1, . . . , 6, are the solution of the system. If not all of the conditions in Step 13 are verified, then set k = k + 1 and return to Step 2. 4.5 Preliminary results and discussion Numerical simulations are run to test the model outcomes in the control state, which corresponds to the condition of a healthy individual. At this stage, the model does not include any autoregulation mechanism and the level of retinal blood flow is fixed to the healthy value of Q = 6.08 · 10−4 mL/s obtained from the coupled model (Chapter 3). The model predictions of partial pressure of oxygen in the blood pc , blood oxygen saturation SO2 in the six compartments of the microvasculature, and partial pressure of oxygen in the retinal tissue Layers 1 − 6 are depicted in Figure 4.5, Figure 4.6, and Figure 4.7, respectively. The qualitative distribution of blood oxygen saturation is analogous to the results in [45, 131] that show a steep decrease in saturation in the capillary compartment. The value of venous oxygen saturation predicted by the model is 61% and is consistent with the value of 60% ± 5.7% measured in [32] (Figure 4.5). The partial pressure of oxygen and the oxygen saturation in the two capillary compartments show a slightly different behavior (Figure 4.5b and 4.6b) due to the diverse level of oxygen partial pressure in the tissue Layers 1 and 3 (Figure 4.7). The model predictions of tissue oxygen partial pressure (Figure 4.7) are in agreement with previous model results in [127]. The improved vascular architecture and oxygen exchange condition in the model allows a heterogeneous distribution of oxygen in the two capillary compartments as 103 (b) blood oxygen partial pressure [mmHg] blood oxygen partial pressure [mmHg] (a) 110 SA LA 100 90 80 70 60 C3 50 C1 40 SV LV 30 20 0 0.5 1 1.5 2 110 100 90 80 70 60 50 C3 40 C1 30 20 1.02 distance [cm] 1.03 1.04 1.05 1.06 distance [cm] Figure 4.5. (a) Model predictions of blood oxygen partial pressure as a function of position in the network in the six compartments LA,SA,C1,C3,SV,LV at control state.(b) Blood oxygen partial pressure as a function of position in the network in the C1 and C3 compartments. shown in Figure 4.5b and 4.6b, dictated by the local oxygen source and consumption in the different retinal layers. Furthermore the improved retinal anatomy adds additional levels of detail to the descriptive ability of the model, allowing for separate contributions of the inner and outer retina. Even though the results presented here are preliminary, they stress the importance of the improvements introduced by the retinal anatomy, vascular architecture and oxygen exchange in the model. Thus, the coupling between this model, the autoregulatory mechanisms introduced in Chapter 3, and the macrocirculation model introduced in Chapter 2 could help to interpret the outcomes of clinical studies and lead to a better understanding of ocular diseases. 104 (a) (b) 1 1 blood oxygen saturation SA 0.95 blood oxygen saturation LA 0.95 0.9 0.85 0.8 0.75 C3 0.7 C1 0.65 SV 0.9 0.85 0.8 C3 0.75 C1 0.7 0.65 LV 0.6 0.6 0.55 0 0.5 1 1.5 0.55 2 1.02 distance [cm] 1.03 1.04 1.05 1.06 distance [cm] Figure 4.6. (a) Model predictions of blood oxygen saturation as a function of position in the network in the six compartments LA,SA,C1,C3,SV,LV at control state. (b) Blood oxygen saturation as a function of position in the network in the C1 and C3 compartments. 80 2 tissue oxygen partial pressure [mmHg] 1 4 3 5 6 70 60 50 40 30 20 10 0 0 50 100 150 200 250 retinal depth [µm] Figure 4.7. Model predictions of tissue oxygen partial pressure distribution in the six retinal layers in Figure 4.1. 105 5. MATHEMATICAL MODELING OF AQUEOUS HUMOR FLOW IOP is an important input parameter in the models presented in Chapters 2 and 3. In those models, the value of IOP is set a priori and is varied to define different clinical states. These models did not include the various factors within the eye that combine to yield a specific IOP value. In this Chapter, we move a first step towards the study of the multiple mechanisms that contribute to determine the IOP level. The flow of AH inside the eye plays an important role in determining the level of IOP [137–139]. Thus, we compute IOP as the solution of a simplified mathematical model describing the balance between AH production and drainage; we then perform a sensitivity analysis to quantify the influence of variations in physiological parameters on the IOP level in different situations of clinical interest. The model accounts for the variability of these parameters in a systematic manner in order to identify patientspecific factors that may influence the efficacy of IOP-lowering medications, with the ultimate goal of contributing to the development of novel, effective, and individualized therapeutic approaches in glaucoma management. In Section 5.1 we describe the model and its governing equations. In Section 5.2 we present the model results for four different clinical cases of medical interest. In Section 5.3 we discuss and interpret the model outcomes, comparing them with the results of different clinical studies. The material and results presented in this Chapter are partially published in [140]. 5.1 Methods To analyze AH flow, we implement a mathematical model that describes the steady state value of IOP as the result of the balance between AH production and drainage. 106 cBP Gin Guv IOP EVP Gtm J Figure 5.1. Network model for production and drainage of aqueous humor. Changes in ocular blood volume, mainly localized in the choroid, are conjectured to affect the time variations of IOP [138], but they are not considered here. The flow Q of AH is described using the analogy between the flow of a fluid in a hydraulic network and the flow of a current in an electric circuit, depicted in Figure 5.1. Variable conductances are indicated with an arrow. The flow of AH Q through a conductor is described using Ohm’s law: Q = G∆P, (5.1) where the conductance G (= 1/R, reciprocal of resistance) represents the ability to allow the passage of fluid through a physiological pathway, and ∆P represents the pressure drop along that pathway. The contribution of protein and low-molecular components to AH production is included in the model through a current source of intensity equal to J, which thus represents an AH flow rate. AH is produced at the level of the ciliary body by a combination of a passive mechanism (the ultrafiltration) and an active mechanism (the ionic secretion) and is modulated by the total inflow facility Gin [137, 139, 141]. Here, the term facility 107 indicates hydraulic conductance, namely a flow rate per units of pressure. The total flow Jin of AH entering the eye is given by Jin = Juf + Jsecr , (5.2) where Juf and Jsecr are the flows due to ultrafiltration and active secretion, respectively. The ultrafiltration from the ciliary circulation consists of flow of transparent fluid across semipermeable membranes (including vascular walls, stroma and epithelial cells) and is driven by differences in hydrostatic pressures (cBP-IOP), where cBP is the blood pressure in the capillaries of the ciliary body, and oncotic pressures (∆Πp ) modulated by a protein reflection coefficient (σp ). We thus write Juf as   Juf = Gin (cBP − IOP) − σp ∆Πp . (5.3) The inflow due to active ionic secretion is proportional to the blood/AH osmotic pressure difference (∆Πs ) via a reflection coefficient for low-molecular components (σs ), and it is similarly written as Jsecr = Gin [−σs ∆Πs ] . (5.4) Thus, the value of AH flow generated by J (Figure 5.1) is equal to  J = −Gin σp ∆Πp + Jsecr = −Gin σp ∆Πp + σs ∆Πs . (5.5) The drainage of AH from the eye is driven by passive mechanisms through two different pathways. The trabecular pathway, also known as the conventional pathway, consists of AH flow through the trabecular meshwork into the Schlemm’s canal and the episcleral veins. The uveoscleral pathway, also known as the non-conventional pathway, consists of AH flow through the ciliary muscle and into the supraciliary space. Thus, the total flow Jout of AH leaving the eye is given by Jout = Jtm + Juv , (5.6) where Jtm and Juv are the flows via the trabecular and uveoscleral pathways, respectively. As proposed by Brubaker [142], the trabecular pathway model consists 108 of flow through a nonlinear resistor positioned between the anterior chamber (where the pressure is equal to IOP) and the episcleral veins (where the pressure is equal to the episcleral venous pressure (EVP)) with outflow facility Gtm , and is given by the following equation: Jtm = Gtm (IOP − EVP), with Gtm = 1  , R0 1 + Θ (IOP − EVP) (5.7) where R0 is the resistance when IOP equals EVP, and Θ is the outflow obstruction coefficient. The contribution of the uveoscleral pathway is modeled as the flow through a nonlinear resistor connected to the ground [139], with outflow facility Guv depending nonlinearly on the pressure through a Michaelis-Menten type relation [143]: Juv = Guv (IOP − 0), with Guv = k1 , k2 + IOP (5.8) where k1 is the maximum value attainable by the uveoscleral flow rate and k2 is the Michaelis constant for the uveoscleral flow rate, namely the pressure value for which the uveoscleral flow rate is half of k1 . The value of the pressure at the internal node of the system depicted in Figure 5.1 corresponds to the value of IOP. The steady state value of IOP, resulting from the balance between production and drainage of AH, namely Jin = Jout , can be written as: Juf + Jsecr = Jtm + Juv , (5.9) or, equivalently:   Gin (cBP − IOP) − σp ∆Πp − σs ∆Πs = 1 k1   (IOP − EVP) + = IOP. k2 + IOP R0 1 + Θ (IOP − EVP) (5.10) This is a scalar third-order polynomial equation in one unknown, IOP, that can be found explicitly. Control state values for the parameters, defined to represent typical conditions of a healthy eye, are indicated with an overline bar in Table 5.1. To include potential sources of uncertainties as well as to identify and rank parameters having the most important influence on IOP, we applied a global stochastic 109 Table 5.1. Control state values for the parameters in the model for aqueous humor flow in equation (5.10). Parameter µL min mmHg Blood pressure in the capillaries Total inflow facility Gin , Value Source 0.3 [141] 27.5 [138, 139, 141] of the ciliary body cBP, mmHg Blood/AH oncotic pressure difference ∆Πp , mmHg 25 [141] Reflection coefficient for proteins σp 1 [141] Blood/AH osmotic pressure difference ∆Πs , mmHg −450 [141] Reflection coefficient for low-molecular components σs 0.0515 [141] Episcleral venous pressure EVP, mmHg 8 [139] 2.2 [142] Trabecular outflow obstruction coefficient Θ, mmHg−1 0.012 [142] Maximum uveoscleral flow rate k1 , µL/min 0.4 [139] 5 [139] Trabecular outflow resistance (when pressure gradient equals 0) R0 , mmHg min/µL Pressure at which uveoscleral flow rate is at half maximum k2 , mmHg sensitivity analysis to the model described above. We considered stochastic variations in cBP following a normal distribution, and in Gin , ∆Πs , G0 = 1/R0 , k1 and EVP following a uniform distribution. By using the probability distribution of IOP, we computed variance-based sensitivity indices, also known as Sobol indices [144], and the probability density function [145], which describes the relative frequency of a given IOP value. In sensitivity analysis, Sobol indices explain the importance of an input factor on the variance of the output (that corresponds to the level of IOP). For each parameter, its direct influence on IOP is quantified in terms of first-order Sobol 110 indices, and the influence through interactions with other parameters is identified by means of the total Sobol indices. The values of first-order and total indices can be estimated via Monte Carlo simulations [144], or via reduced order models using polynomial chaos expansion [146]. The former method is very costly from the computational viewpoint as it requires many evaluations to ensure convergence, whereas the latter requires considerably less evaluations. Both methods have been compared and provide similar results. Here, we show the results obtained using the polynomial chaos reduced model. 5.2 Results The model is used to compute the IOP distribution in four different cases of clini- cal interest: (i) ocular normotensive healthy subjects (ONT); (ii) ocular hypertensive subjects (OHT); (iii) ONT subjects treated with IOP-lowering medications (ONTm); and (iv) OHT subjects treated with IOP-lowering medications (OHTm). The IOP probability density function and first and total Sobol indices are reported in Figure 5.2 for the four cases. Mean values, standard deviations, and skewness of the IOP distribution in the four cases are reported in Table 5.2. Model simulations and results are described below. (i) ONT subjects. The mean values of cBP, Gin , ∆Πs , G0 , k1 and EVP are set equal to their control state values and are summarized in Table 5.1. Variations in cBP are deduced from variations in MAP. Specifically, we write cBP= α MAP, where α = 0.296 is chosen as to obtain cBP= 27.5 mmHg when MAP= 93 mmHg; we assumed a normal distribution for MAP of 93 ± 7.6 mmHg [148]. Variations in Gin , ∆Πs , G0 , k1 and EVP are assumed to follow a uniform distribution with a variation of 15% with respect to their control state values. Simulation outcomes: The IOP probability density function for ONT subjects (Figure 5.2a) fits a right-skewed Gaussian curve with a frequency peak of 25% at 15.13 mmHg and a skewness of 0.2, which is in a very good agreement with the 111 Table 5.2. Mean values, standard deviations and skewness of the distribution of intraocular pressure (IOP) resulting from the sensitivity analysis of the mathematical model in equation (5.10) for four cases of clinical interest. Skewness of IOP Clinical cases IOP, mmHg distribution Ocular normotensive (ONT) 15.13 ± 1.58 0.2 Ocular hypertensive (OHT) 20.12 ± 2.35 0.09 12.58 ± 1.32 0.17 15.81 ± 2.03 0.08 Ocular normotensive treated with IOP-lowering medications (ONTm) Ocular hypertensive treated with IOP-lowering medications (OHTm) results from a population-based study on approximately 12, 000 subjects [147] (green curve in Figure 5.2a). The results for the Sobol indices (Figure 5.2b) suggest that IOP is strongly influenced by cBP and ∆Πs and mildly influenced by the levels of Gin , G0 and EVP. The influence of k1 on IOP appears to be minimal. (ii) OHT subjects. OHT condition is simulated by decreasing the mean value of the trabecular meshwork outflow facility as suggested by several clinical observations [149, 150]. Thus, here we set G0 = 0.3G0 , leaving the mean values of the other parameters at control state values. Simulation outcomes: IOP probability density function in the OHT case (Figure 5.2c) fits a Gaussian curve, but with a frequency peak of 15% at 20.12 mmHg and with a more symmetric profile than the ONT Gaussian curve (skewness = 0.09). The Sobol indices values for OHT subjects (Figure 5.2d) show a stronger dependence of IOP on cBP and ∆Πs and a weaker dependence of IOP on Gin , 112 G0 and EVP than for ONT subjects. The influence of k1 on IOP remains minimal. (iii) ONT subjects treated with IOP-lowering medications (ONTm). We model the effect of IOP-lowering medications by reducing the active ionic secretion by 25%, which sets the mean value of the blood/AH osmotic pressure difference to = 0.75∆Πs ; the mean values of the other parameters remained at control state. This modeling choice is justified by the fact that the sensitivity analyses in both the ONT and OHT cases have identified ∆Πs as an important determinant of IOP levels; in addition, clinical evidence and studies also support this notion [137–139]. Simulation outcomes: The IOP probability density function in the ONTm case (Figure 5.2e) fits a right-skewed Gaussian curve with a frequency peak of 30% at 12.58 mmHg and a skewness of 0.08. Thus, our simulations predict a reduction of 2.55 mmHg in the mean value of IOP when IOP-lowering medications are administered to ONT subjects. The results of the Sobol indices (Figure 5.2f) suggest that IOP is strongly influenced by cBP and ∆Πs and mildly influenced by the levels of Gin , G0 and EVP. The influence of k1 on IOP is again minimal. (iv) OHT subjects treated with IOP-lowering medications (OHTm). We simultaneously account for OHT conditions and IOP-lowering treatment by setting the mean values of G0 and ∆Πs to G0 = 0.3G0 and ∆Πs = 0.75∆Πs , leaving the mean values of the other parameters at control state values. Simulation outcomes: IOP probability density function in the OHTm case (Figure 5.2g) fits a Gaussian curve with a frequency peak of 20% at 15.81 mmHg and has a more symmetric profile than the curve in the ONTm case (skewness = 0.08). Thus, our simulations predict a reduction of 4.31 mmHg in the mean value of IOP when IOP-lowering medications are administered to OHT subjects. The results on Sobol indices (Figure 5.2h) are similar to those obtained in the ONTm case, but with a weaker contribution of Gin , G0 and EVP. 113 Our results demonstrate that first-order and total Sobol indices do not present noticeable differences in any of the four simulated scenarios, suggesting that higher order interactions among the selected factors are minimal. 5.3 Discussion The model reproduced conditions of normal ocular tension, with blood pressure and IOP values within physiological ranges, and was subsequently used to simulate the effect of IOP-lowering medications. The proposed model suggests that the outcomes of IOP-lowering treatments depend on the initial IOP level of the patient and on the patient’s individual clinical condition. Specifically, the model predicts mean IOP reductions of 2.55 mmHg and 4.31 mmHg when the pre-treatment IOP mean values are 15.13 mmHg and 20.12 mmHg, respectively. These predictions are in good agreement with Rulo et al. [151] who reported mean IOP reductions of 2 mmHg and 4.2 mmHg for pre-treatment IOP mean values of 15.3 mmHg and 18.4 mmHg, respectively. However, it is important to remark that the study by Rulo et al. utilized Latanoprost, a prostaglandin analog that increases AH drainage, whereas we modeled IOP-lowering medications by decreasing AH production. Other studies reported IOP reductions ranging from 3 mmHg to 4.4 mmHg in response to brinzolamide [152], from 4.5 mmHg to 6.1 mmHg in response to dorzolamide [153], and from 2.4 mmHg to 4.5 mmHg in response to Latanoprost [154]. The mean IOP reductions reported in these studies [155] are close or slightly higher than those predicted by our model; this might be due to the fact that these studies started from higher pre-treatment IOP levels (ranging from 23.8 mmHg to 28.9 mmHg) than those considered in our simulations. Our analysis also suggests that IOP-lowering effects are more pronounced when AH production is affected rather than AH drainage. The effects of lowering IOP are also more apparent when trabecular outflow is increased instead of the uveoscleral outflow. Another interesting finding of our analysis is that a patient’s blood pres- 114 sure strongly influences the outcomes of IOP-lowering treatments, which may explain why the effect of some drugs differ between day-time and night-time and/or among individuals [151, 154, 156, 157]. A further investigation that incorporates a theoretical model coupling AH production and drainage with ocular blood flow may lead to a better understanding of this delicate, yet important, relationship [46, 139, 158]. : Clinical data [147] Relative frequency 0.25 0.25 0.20 0.20 Sobol Indices Total Sobol Indices (b): ONT 0.75 0.75 0.50 0.50 0.15 0.15 0.10 0.10 0.25 0.25 0.05 0.05 0.00 0.00 5 5 10 10 15 20 25 15 p_ant20 [mmHg] 25 IOP [mmHg] 0.35 0.35 30 30 35 35 0.00 0.00 1.00 1.00 (c): OHT 0.30 0.30 MAP Lin BP Gin r0 G0 k1 EVP ∆πs k1 EVP ∆Πs Sobol Indices Total Sobol Indices (d): OHT 0.75 0.75 0.25 0.25 Relative frequency IOP relative frequency 1.00 1.00 (a): ONT 0.30 0.30 0.20 0.20 0.50 0.50 0.15 0.15 0.10 0.10 0.25 0.25 0.05 0.05 0.00 0.00 5 5 10 10 15 20 25 15 p_ant20 [mmHg] 25 IOP [mmHg] 0.35 0.35 30 30 35 35 (e): ONTm 0.30 0.30 0.00 0.00 1.00 1.00 MAP Lin BP Gin r0 G0 k1 EVP ∆πs k1 EVP ∆Πs Sobol Indices Total Sobol Indices (f): ONTm 0.75 0.75 0.25 0.25 Relative frequency IOP relative frequency IOP relative frequency 0.35 0.35 0.20 0.20 0.50 0.50 0.15 0.15 0.10 0.10 0.25 0.25 0.05 0.05 0.00 0.00 5 5 10 10 15 20 25 15 p_ant20 [mmHg] 25 IOP [mmHg] 0.35 0.35 30 30 35 35 (g): OHTm 0.30 0.30 0.00 0.00 1.00 1.00 MAP Lin BP Gin r0 G0 k1 EVP ∆πs k1 EVP ∆Πs Sobol Indices Total Sobol Indices (h): OHTm 0.75 0.75 0.25 0.25 Relative frequency IOP relative frequency 115 0.20 0.20 0.50 0.50 0.15 0.15 0.10 0.10 0.25 0.25 0.05 0.05 0.00 0.00 5 5 10 10 15 20 25 15 p_ant20 [mmHg] 25 IOP [mmHg] 30 30 35 35 0.00 0.00 MAP Lin BP Gin r0 G0 k1 EVP ∆πs k1 EVP ∆Πs Figure 5.2. Probability density function of intraocular pressure (IOP) (left figures), and Sobol indices resulting from the sensitivity analysis performed on the mathematical model of equation (5.10) when variations in ciliary capillary blood pressure cBP, total inflow facility Gin , blood/AH osmotic pressure difference ∆Πs , trabecular outflow facility G0 , uveoscleral outflow facility k1 and episcleral venous pressure EVP are considered (right figures). The model predictions are obtained for the clinical cases ONT, OHT, ONTm, OHTm. 116 6. CONCLUSIONS In this thesis, mathematical models for the retinal circulation, the transport of oxygen in the retina, and the production and drainage of AH are proposed. Summarizing: • We developed a compartmental time-dependent mathematical model for the retinal circulation that accounts for blood flow in the central retinal vessels and microvasculature, vascular regulation, and the effect of IOP on the retinal vasculature. The model has been used to compute the values of retinal blood flow, blood velocity and pressure distribution in the retinal vasculature for six theoretical patients with different clinical conditions. The model results have been compared with clinical data to clarify the impact of changes in IOP, BP and AR on retinal hemodynamics. The model has been directly applied to the interpretation of the outcomes of trabeculectomy (Chapter 2). • We used mathematical methods to describe the difference between compressible and collapsible tubes and implemented these methods to describe the passive resistances of the CRA, venules and CRV compartments (Section 2.1.5). • We proposed a coupling between the model described in Chapter 2 and the mechanistic model for autoregulation in [45]. The coupling introduced the Krogh-cylinder type model for the diffusion of oxygen from the retinal vasculature into the retinal tissue and four vascular regulation mechanisms that substituted the phenomenological approach used in Chapter 2. This improved model has been used to predict the level of tissue oxygenation for alterations in BP and tissue oxygen demand. Furthermore the model predictions revealed the relative importance of the different regulation mechanisms to achieve VR (Chapter 3). 117 • We improved the description of the transport of oxygen in the retinal microvasculature and retinal tissue by including a more physiological representation of retinal anatomy, vascular architecture and oxygen exchange in Chapter 4. The preliminary results obtained with this model for oxygen blood saturation in the retinal microvasculature and partial pressure of oxygen in the retinal tissue are in good agreement with the literature. • We implemented a mathematical model that describes the production and drainage of AH, which affect the level of IOP. A sensitivity analysis was preformed to quantify the effect of some physiological parameters on the IOP distribution in theoretical patients of clinical interest and to understand the uncertainty of the outcomes of clinical studies on the effect of IOP-lowering medications. The interpretation of the model results contributes to answer the following questions of clinical relevance: How do IOP variations affect retinal hemodynamics? The model predictions in Chapter 2 show that the hemodynamic response of the system to changes in IOP varies depending on the level of BP and the functionality of AR. The AR plateau is influenced by the level of BP, occurring for higher values of pressure for hypertensive individuals. Thus, different values of retinal blood flow can correspond to the same value of IOP depending on the level of BP and functionality of AR. The model results for CRA blood velocity suggest that the same reduction in IOP can lead to small or large changes in velocity depending on the initial value of IOP. In addition, the model predictions show that, for some theoretical patients, the IOP-interval for which the retinal blood flow is almost constant does not overlap with the IOP-interval where PSV and EDV are almost constant. This results is interesting because, in clinical studies, the measured value of velocity is used as a surrogate to 118 compute the level of blood flow. However, the model results show that in some cases this will not yield an accurate estimate of the trends in retinal blood flow. The model predictions for trabeculectomy give a similar picture, in which the outcomes of a decrease in IOP strongly depend on the level of BP and AR functionality of the patient. What is the effect of collapsible tubes? The collapsibility of veins plays a significant role on the pressure distribution in the system. As IOP increases, the CRV and the venules compartments collapse, causing an increase in the upstream pressure. This passive feedback mechanism, consistent with experimental observations, helps the retinal microvasculature to withstand the increased IOP and improves the level of perfusion. What is the effect of the coupling with [45] and what mechanisms play a major role in vascular regulation? The results in Chapter 3 stress the importance of coupling the macrocirculation model with the microcirculation model. Although the coupled and microcirculation models have a similar behavior in the range of pressures where VR works properly, the results notably differ outside this range. This suggests that the coupled model provides a more comprehensive and accurate depiction of hemodynamics in the retina than the sole microcirculation model. When comparing the contribution of the vascular regulation mechanisms, the model highlights the importance of the metabolic and CO2 responses. Impairing either of these mechanisms could cause insufficient tissue oxygenation and therefore ischemic damage. The model also predicts that blood flow will increase more for an increase in tissue oxygen demand than for the same percentage decrease in O2 demand. How does retinal anatomy, vascular architecture and oxygen exchange affect tissue oxygenation? 119 The model derived in Chapter 4 includes a heterogeneous description of the retinal tissue, adapts the vascular network to accommodate for two distinct levels of capillaries, and introduces a new model for the exchange of oxygen between retinal tissue and vasculature that responds to local changes in oxygen partial pressure. This more realistic tissue/vascular architecture, together with the improved model for oxygen exchange, isolate the contributions of the inner and outer retina in the oxygen profiles predicted by the model, enhancing the model’s predictive ability. What are the major factors influencing the level of IOP? The statistical analysis performed in Chapter 5 suggests that the inconsistencies of the effect of IOP-lowering medication might depend on the level of IOP before the medicine was administrated. Furthermore a greater decrease in IOP is predicted when the medicine affects AH production rather then AH drainage. The statistical analysis reveals that BP is one of the primary factors that influences the value of IOP, which might explain the difference in effectiveness of medications among individuals. What insight can mathematical modeling provide to clinicians? The results of this thesis demonstrate that mathematical modeling can be used alongside clinical and animal studies as (i) a tool to understand the complex relationship among different factors at play, (ii) a virtual lab where theories and treatments can be tested on theoretical patients and guide clinical developments, and (iii) an instrument to quantify or estimate factors that are difficult or impossible to measure. REFERENCES 120 REFERENCES [1] H. Gray, Anatomy of the human body. Lea & Febiger, 1918. [2] A. Harris, L. Kagemann, R. Ehrlich, C. Rospigliosi, D. Moore, and B. Siesky, “Measuring and interpreting ocular blood flow and metabolism in glaucoma,” Canadian Journal of Ophthalmology/Journal Canadien d’Ophtalmologie, vol. 43, no. 3, pp. 328–336, 2008. [3] M. C. Leske, “Open-angle glaucoma—an epidemiologic overview,” Ophthalmic epidemiology, vol. 14, no. 4, pp. 166–172, 2007. [4] M. C. Leske, A. Heijl, L. Hyman, B. Bengtsson, L. Dong, Z. Yang, E. Group et al., “Predictors of long-term progression in the early manifest glaucoma trial,” Ophthalmology, vol. 114, no. 11, pp. 1965–1972, 2007. [5] R. Ehrlich, A. Harris, N. S. Kheradiya, D. M. Winston, T. A. Ciulla, and B. Wirostko, “Age-related macular degeneration and the aging eye,” Clinical interventions in aging, vol. 2008, p. 473, 2008. [6] B. Pemp and L. Schmetterer, “Ocular blood flow in diabetes and age-related macular degeneration,” Canadian Journal of Ophthalmology/Journal Canadien d’Ophtalmologie, vol. 43, no. 3, pp. 295–301, 2008. [7] S. Rassam, V. Patel, and E. Kohner, “The effect of experimental hypertension on retinal vascular autoregulation in humans: a mechanism for the progression of diabetic retinopathy,” Experimental physiology, vol. 80, no. 1, pp. 53–68, 1995. [8] G. Alagöz, K. Gürel, A. Bayer, D. Serin, S. Çelebi, and Ş. Kükner, “A comparative study of bimatoprost and travoprost: effect on intraocular pressure and ocular circulation in newly diagnosed glaucoma patients,” Ophthalmologica, vol. 222, no. 2, pp. 88–95, 2008. [9] O. Findl, K. Strenn, M. Wolzt, R. Menapace, C. Vass, H.-G. Eichler, and L. Schmetterer, “Effects of changes in intraocular pressure on human ocular haemodynamics,” Current eye research, vol. 16, no. 10, pp. 1024–1029, 1997. [10] F. Galassi, B. Giambene, A. Corvi, G. Falaschi, and U. Menchini, “Retrobulbar hemodynamics and corneal surface temperature in glaucoma surgery,” International ophthalmology, vol. 28, no. 6, pp. 399–405, 2008. [11] A. Harris, K. Joos, M. Kay, D. Evans, R. Shetty, W. E. Sponsel, and B. Martin, “Acute iop elevation with scleral suction: effects on retrobulbar haemodynamics.” British journal of ophthalmology, vol. 80, no. 12, pp. 1055–1059, 1996. 121 [12] K. K. Huber-van der Velden, A. Lux, K. Severing, M. K. J. Klamann, S. Winterhalter, and A. Remky, “Retrobulbar hemodynamics before and after oculopression with and without dorzolamide,” Current eye research, vol. 37, no. 8, pp. 719–725, 2012. [13] O. G. Koz, A. Ozsoy, A. Yarangumeli, S. K. Kose, and G. Kural, “Comparison of the effects of travoprost, latanoprost and bimatoprost on ocular circulation: a 6-month clinical trial,” Acta Ophthalmologica Scandinavica, vol. 85, no. 8, pp. 838–843, 2007. [14] A. Synder, E. Augustyniak, I. Laudańska-Olszewska, and A. Wesołek-Czernik, “[evaluation of blood-flow parameters in extraocular arteries in patients with primary open-angle glaucoma before and after surgical treatment].” Klinika oczna, vol. 106, no. 1-2 Suppl, pp. 206–208, 2003. [15] J. R. Trible, R. C. Sergott, G. L. Spaeth, R. P. Wilson, L. J. Katz, M. R. Moster, and C. M. Schmidt, “Trabeculectomy is associated with retrobulbar hemodynamic changes: a color doppler analysis,” Ophthalmology, vol. 101, no. 2, pp. 340–351, 1994. [16] C. Akarsu, S. Yılmaz, P. Taner, and A. Ergin, “Effect of bimatoprost on ocular circulation in patients with open-angle glaucoma or ocular hypertension,” Graefe’s Archive for Clinical and Experimental Ophthalmology, vol. 242, no. 10, pp. 814–818, 2004. [17] O. K. Arikan, C. Akarsu, B. Unal, A. Ergin, and C. Koç, “Effect of oxymetazoline nasal spray on intraocular pressure and retrobulbar hemodynamics.” Journal of otolaryngology, vol. 35, no. 1, 2006. [18] L. B. Cantor, “The effect of trabeculectomy on ocular hemodynamics.” Transactions of the American Ophthalmological Society, vol. 99, p. 241, 2001. [19] M. L. Conway, M. Wevill, A. Benavente-Perez, and S. L. Hosking, “Ocular blood-flow hemodynamics before and after application of a laser in situ keratomileusis ring,” Journal of Cataract & Refractive Surgery, vol. 36, no. 2, pp. 268–272, 2010. [20] E. Erkin, S. Tarhan, O. Kayikçioğlu, H. Deveci, C. Güler, and C. Göktan, “Effects of betaxolol and latanoprost on ocular blood flow and visual fields in patients with primary open-angle glaucoma.” European journal of ophthalmology, vol. 14, no. 3, pp. 211–219, 2003. [21] G. Fuchsjäger-Mayrl, M. Georgopoulos, A. Hommer, G. Weigert, B. Pemp, C. Vass, G. Garhöfer, and L. Schmetterer, “Effect of dorzolamide and timolol on ocular pressure: blood flow relationship in patients with primary open-angle glaucoma and ocular hypertension,” Investigative ophthalmology & visual science, vol. 51, no. 3, pp. 1289–1296, 2010. [22] A. Hommer, P. Sperl, H. Resch, A. Popa-Cherecheanu, C. Qiao, L. Schmetterer, and G. Garhöfer, “A double-masked randomized crossover study comparing the effect of latanoprost/timolol and brimonidine/timolol fixed combination on intraocular pressure and ocular blood flow in patients with primary open-angle glaucoma or ocular hypertension,” Journal of Ocular Pharmacology and Therapeutics, vol. 28, no. 6, pp. 569–575, 2012. 122 [23] A. Harris, H. J. Garzozi, L. McCranor, E. Rechtman, C.-W. Yung, and B. Siesky, “The effect of latanoprost on ocular blood flow,” International ophthalmology, vol. 29, no. 1, pp. 19–26, 2009. [24] Ü. Ü. Inan, S. S. Ermis, A. Yücel, and F. Öztürk, “The effects of latanoprost and brimonidine on blood flow velocity of the retrobulbar vessels: a 3-month clinical trial,” Acta Ophthalmologica Scandinavica, vol. 81, no. 2, pp. 155–160, 2003. [25] C. James, “Effect of trabeculectomy on pulsatile ocular blood flow.” British journal of ophthalmology, vol. 78, no. 11, pp. 818–822, 1994. [26] M. Kaup, N. Plange, M. Niegel, A. Remky, and O. Arend, “Effects of brinzolamide on ocular haemodynamics in healthy volunteers,” British journal of ophthalmology, vol. 88, no. 2, pp. 257–262, 2004. [27] D. Poinoosawmy, A. Indar, C. Bunce, D. Garway-Heath, and R. Hitchings, “Effect of treatment by medicine or surgery on intraocular pressure and pulsatile ocular blood flow in normal-pressure glaucoma,” Graefe’s Archive for Clinical and Experimental Ophthalmology, vol. 240, no. 9, pp. 721–726, 2002. [28] A. Stankiewicz, M. Misiuk-Hojło, I. Grabska-Liberek, B. Romanowska-Dixon, J. Wierzbowska, J. Wasyluk, M. Mulak, I. Szuścik, J. Sierdziński, R. Ehrlich et al., “Intraocular pressure and ocular hemodynamics in patients with primary open-angle glaucoma treated with the combination of morning dosing of bimatoprost and dorzolamide hydrochloride,” Acta ophthalmologica, vol. 89, no. 1, pp. e57–e63, 2011. [29] G. T. Dorner, G. Garhöfer, K. H. Huemer, C. E. Riva, M. Wolzt, and L. Schmetterer, “Hyperglycemia affects flicker-induced vasodilation in the retina of healthy subjects,” Vision Research, vol. 43, no. 13, pp. 1495–1500, 2003. [30] G. Garhöfer, C. Zawinka, H. Resch, K. Huemer, G. Dorner, and L. Schmetterer, “Diffuse luminance flicker increases blood flow in major retinal arteries and veins,” Vision Research, vol. 44, no. 8, pp. 833–838, 2004. [31] G. Garhöfer, C. Zawinka, H. Resch, P. Kothy, L. Schmetterer, and G. Dorner, “Reduced response of retinal vessel diameters to flicker stimulation in patients with diabetes,” British Journal of Ophthalmology, vol. 88, no. 7, pp. 887–891, 2004. [32] M. Hammer, W. Vilser, T. Riemer, F. Liemt, S. Jentsch, J. Dawczynski, and D. Schweitzer, “Retinal venous oxygen saturation increases by flicker light stimulation,” Investigative Ophthalmology & Visual Science, vol. 52, no. 1, pp. 274– 277, 2011. [33] A. Mandecka, J. Dawczynski, M. Blum, N. Müller, C. Kloos, G. Wolf, W. Vilser, H. Hoyer, and U. A. Müller, “Influence of flickering light on the retinal vessels in diabetic patients,” Diabetes Care, vol. 30, no. 12, pp. 3048–3052, 2007. [34] S. H. Hardarson, S. Basit, T. E. Jonsdottir, T. Eysteinsson, G. H. Halldorsson, R. A. Karlsson, J. M. Beach, J. A. Benediktsson, and E. Stefansson, “Oxygen saturation in human retinal vessels is higher in dark than in light,” Investigative Ophthalmology & Visual Science, vol. 50, no. 5, pp. 2308–2311, 2009. 123 [35] Y.-Y. I. Shih, L. Wang, B. H. De La Garza, G. Li, G. Cull, J. W. Kiel, and T. Q. Duong, “Quantitative retinal and choroidal blood flow during light, dark adaptation and flicker light stimulation in rats using fluorescent microspheres,” Current Eye Research, vol. 38, no. 2, pp. 292–298, 2013. [36] N. D. Wangsa-Wirawan and R. A. Linsenmeier, “Retinal oxygen: fundamental and clinical aspects,” Archives of Ophthalmology, vol. 121, no. 4, pp. 547–557, 2003. [37] J. Caprioli, A. L. Coleman et al., “Blood pressure, perfusion pressure, and glaucoma,” American journal of ophthalmology, vol. 149, no. 5, pp. 704–712, 2010. [38] S. S. Hayreh, M. B. Zimmerman, P. Podhajsky, and W. L. Alward, “Nocturnal arterial hypotension and its role in optic nerve head and ocular ischemic disorders,” American journal of ophthalmology, vol. 117, no. 5, pp. 603–624, 1994. [39] Z. He, C. T. Nguyen, J. A. Armitage, A. J. Vingrys, and B. V. Bui, “Blood pressure modifies retinal susceptibility to intraocular pressure elevation,” PLoS One, vol. 7, no. 2, p. e31104, 2012. [40] D. Leighton and C. Phillips, “Systemic blood pressure in open-angle glaucoma, low tension glaucoma, and the normal eye.” The British journal of ophthalmology, vol. 56, no. 6, p. 447, 1972. [41] F. Memarzadeh, M. Ying-Lai, J. Chung, S. P. Azen, and R. Varma, “Blood pressure, perfusion pressure, and open-angle glaucoma: the los angeles latino eye study,” Investigative ophthalmology & visual science, vol. 51, no. 6, pp. 2872–2877, 2010. [42] F. Galassi, B. Giambene, and R. Varriale, “Systemic vascular dysregulation and retrobulbar hemodynamics in normal-tension glaucoma,” Investigative ophthalmology & visual science, vol. 52, no. 7, pp. 4467–4471, 2011. [43] D. Moore, A. Harris, D. WuDunn, N. Kheradiya, and B. Siesky, “Dysfunctional regulation of ocular blood flow: A risk factor for glaucoma?” Clinical ophthalmology (Auckland, NZ), vol. 2, no. 4, p. 849, 2008. [44] D. Sines, A. Harris, B. Siesky, I. Januleviciene, C. Haine, C. Yung, Y. Catoira, and H. Garzozi, “The response of retrobulbar vasculature to hypercapnia in primary open-angle glaucoma and ocular hypertension,” Ophthalmic research, vol. 39, no. 2, pp. 76–80, 2007. [45] J. Arciero, A. Harris, B. Siesky, A. Amireskandari, V. Gershuny, A. Pickrell, and G. Guidoboni, “Theoretical analysis of vascular regulatory mechanisms contributing to retinal blood flow autoregulation,” Investigative ophthalmology & visual science, vol. 54, no. 8, pp. 5584–5593, 2013. [46] G. Guidoboni, A. Harris, S. Cassani, J. Arciero, B. Siesky, A. Amireskandari, L. Tobe, P. Egan, I. Januleviciene, and J. Park, “Intraocular pressure, blood pressure, and retinal blood flow autoregulation: A mathematical model to clarify their relationship and clinical relevance,” Investigative ophthalmology & visual science, vol. 55, no. 7, pp. 4105–4118, 2014. 124 [47] S. Cassani, G. Guidoboni, I. Januleviciene, L. Carichino, B. Siesky, L. A. Tobe, A. Amireskandari, D.-P. Baikstiene, and A. Harris, “3. effect of trabeculectomy on retinal hemodynamics: Mathematical modeling of clinical data,” Integrated Multidisciplinary Approaches in the Study and Care of the Human Eye, p. 29, 2014. [48] G. Guidoboni, A. Harris, L. Carichino, Y. Arieli, and B. A. Siesky, “Effect of intraocular pressure on the hemodynamics of the central retinal artery: a mathematical model.” Mathematical biosciences and engineering: MBE, vol. 11, no. 3, pp. 523–546, 2014. [49] A. Harris, O. Arend, K. Bohnke, E. Kroepfl, R. Danis, and B. Martin, “Retinal blood flow during dynamic exercise,” Graefe’s archive for clinical and experimental ophthalmology, vol. 234, no. 7, pp. 440–444, 1996. [50] M. Iester, P. G. Torre, G. Bricola, A. Bagnis, and G. Calabria, “Retinal blood flow autoregulation after dynamic exercise in healthy young subjects,” Ophthalmologica, vol. 221, no. 3, pp. 180–185, 2007. [51] J. Németh, K. Knézy, B. Tapasztó, R. Kovács, and Z. Harkányi, “Different autoregulation response to dynamic exercise in ophthalmic and central retinal arteries: a color doppler study in healthy subjects,” Graefe’s archive for clinical and experimental ophthalmology, vol. 240, no. 10, pp. 835–840, 2002. [52] W. H. Morgan, D.-Y. Yu, V. A. Alder, S. J. Cringle, R. L. Cooper, P. H. House, and I. J. Constable, “The correlation between cerebrospinal fluid pressure and retrolaminar tissue pressure.” Investigative ophthalmology & visual science, vol. 39, no. 8, pp. 1419–1428, 1998. [53] R. Ren, J. B. Jonas, G. Tian, Y. Zhen, K. Ma, S. Li, H. Wang, B. Li, X. Zhang, and N. Wang, “Cerebrospinal fluid pressure in glaucoma: a prospective study,” Ophthalmology, vol. 117, no. 2, pp. 259–266, 2010. [54] G. T. Dorner, E. Polska, G. Garhöfer, C. Zawinka, B. Frank, and L. Schmetterer, “Calculation of the diameter of the central retinal artery from noninvasive measurements in humans,” Current eye research, vol. 25, no. 6, pp. 341–345, 2002. [55] M. J. Dumskyj, J. E. Eriksen, C. J. Doré, and E. M. Kohner, “Autoregulation in the human retinal circulation: assessment using isometric exercise, laser doppler velocimetry, and computer-assisted image analysis,” Microvascular research, vol. 51, no. 3, pp. 378–392, 1996. [56] G. T. Feke, R. Hazin, C. L. Grosskreutz, and L. R. Pasquale, “Effect of brimonidine on retinal blood flow autoregulation in primary open-angle glaucoma,” Journal of Ocular Pharmacology and Therapeutics, vol. 27, no. 4, pp. 347–352, 2011. [57] G. T. Feke and L. R. Pasquale, “Retinal blood flow response to posture change in glaucoma patients compared with healthy subjects,” Ophthalmology, vol. 115, no. 2, pp. 246–252, 2008. [58] J. E. Grunwald, J. DuPont, and C. E. Riva, “Retinal haemodynamics in patients with early diabetes mellitus.” British journal of ophthalmology, vol. 80, no. 4, pp. 327–331, 1996. 125 [59] T. Takahashi, T. Nagaoka, H. Yanagida, T. Saitoh, A. Kamiya, T. Hein, L. Kuo, and A. Yoshida, “A mathematical model for the distribution of hemodynamic parameters in the human retinal microvascular network,” Journal of biorheology, vol. 23, no. 2, pp. 77–86, 2009. [60] A. Harris, C. Jonescu-Cuypers, L. Kagemann, T. Ciulla, and G. Krieglstein, Atlas of ocular blood flow. Butterworth-Heinemann, 2003. [61] E. J. Lee, T.-W. Kim, and R. N. Weinreb, “Variation of lamina cribrosa depth following trabeculectomylamina cribrosa depth variation,” Investigative ophthalmology & visual science, vol. 54, no. 8, pp. 5392–5399, 2013. [62] J. B. Jonas and L. Holbach, “Central corneal thickness and thickness of the lamina cribrosa in human eyes,” Investigative ophthalmology & visual science, vol. 46, no. 4, pp. 1275–1279, 2005. [63] R. Ren, N. Wang, B. Li, L. Li, F. Gao, X. Xu, and J. B. Jonas, “Lamina cribrosa and peripapillary sclera histomorphometry in normal and advanced glaucomatous chinese eyes with various axial length,” Investigative ophthalmology & visual science, vol. 50, no. 5, pp. 2175–2184, 2009. [64] Y.-c. Fung, Biomechanics: Circulation. 1997. Springer Science & Business Media, [65] A. Quarteroni, M. Tuveri, and A. Veneziani, “Computational vascular fluid dynamics: problems, models and methods,” Computing and Visualization in Science, vol. 2, no. 4, pp. 163–197, 2000. [66] R. L. Armentano, J. G. Barra, J. Levenson, A. Simon, and R. H. Pichel, “Arterial wall mechanics in conscious dogs assessment of viscous, inertial, and elastic moduli to characterize aortic wall behavior,” Circulation Research, vol. 76, no. 3, pp. 468–478, 1995. [67] Y. Fung, “Biomechanics:{M} echanical {P} roperties of {L} iving {T} issues,” 1993. [68] X. Deng and R. Guidoin, “Arteries, veins and lymphatic vessels,” in Handbook of Biomaterial Properties. Springer, 1998, pp. 81–105. [69] D. Baleanu, M. Ritt, J. Harazny, J. Heckmann, R. E. Schmieder, and G. Michelson, “Wall-to-lumen ratio of retinal arterioles and arteriole-to-venule ratio of retinal vessels in patients with cerebrovascular damage,” Investigative ophthalmology & visual science, vol. 50, no. 9, pp. 4351–4359, 2009. [70] E. Wetterer, R. Bauer, and T. Pasch, “Arteriensystem,” in Lehrbuch der Physiologie in Einzeldarstellungen. Springer, 1971, pp. 1–66. [71] A. B. Friedland, “A mathematical model of transmural transport of oxygen to the retina,” Bulletin of mathematical biology, vol. 40, no. 6, pp. 823–837, 1978. [72] H. K. Walker, W. D. Hall, and J. W. Hurst, “Clinical methods,” 1990. [73] A. Strauss and A. Kedra, “Experiences with a new procedure for the measurement of the ophthalmic artery pressure: ophthalmomanometry-doppler.” Medical instrumentation, vol. 21, no. 5, pp. 255–261, 1987. 126 [74] W. D. Lakin, S. A. Stevens, B. I. Tranmer, and P. L. Penar, “A whole-body mathematical model for intracranial pressure dynamics,” Journal of mathematical biology, vol. 46, no. 4, pp. 347–383, 2003. [75] G. Guidoboni, A. Harris, J. C. Arciero, B. A. Siesky, A. Amireskandari, A. L. Gerber, A. H. Huck, N. J. Kim, S. Cassani, and L. Carichino, “Mathematical modeling approaches in the study of glaucoma disparities among people of african and european descents,” Journal of coupled systems and multiscale dynamics, vol. 1, no. 1, p. 1, 2013. [76] J. Jonas, C. Y. Mardin, U. Schlötzer-Schrehardt, and G. Naumann, “Morphometry of the human lamina cribrosa surface.” Investigative ophthalmology & visual science, vol. 32, no. 2, pp. 401–405, 1991. [77] R. E. Norman, J. G. Flanagan, S. M. Rausch, I. A. Sigal, I. Tertinegg, A. Eilaghi, S. Portnoy, J. G. Sled, and C. R. Ethier, “Dimensions of the human sclera: thickness measurement and regional changes with axial length,” Experimental eye research, vol. 90, no. 2, pp. 277–284, 2010. [78] T. Newson and A. El-Sheikh, “Mathematical modeling of the biomechanics of the lamina cribrosa under elevated intraocular pressures,” Journal of biomechanical engineering, vol. 128, no. 4, pp. 496–504, 2006. [79] S.-Y. Woo, A. Kobayashi, W. Schlegel, and C. Lawrence, “Nonlinear material properties of intact cornea and sclera,” Experimental Eye Research, vol. 14, no. 1, pp. 29–39, 1972. [80] A. Quarteroni, L. Formaggia, and N. Ayache, Mathematical modelling and numerical simulation of the cardiovascular system. EPFL, 2002. [81] C. Cancelli and T. Pedley, “A separated-flow model for collapsible-tube oscillations,” Journal of Fluid Mechanics, vol. 157, pp. 375–404, 1985. [82] J. B. Grotberg and O. E. Jensen, “Biofluid mechanics in flexible tubes,” Annual Review of Fluid Mechanics, vol. 36, no. 1, p. 121, 2004. [83] T. J. Pedley, B. S. Brook, and R. S. Seymour, “Blood pressure and flow rate in the giraffe jugular vein,” Philosophical Transactions of the Royal Society of London B: Biological Sciences, vol. 351, no. 1342, pp. 855–866, 1996. [84] A. H. Shapiro, “Steady flow in collapsible tubes,” Journal of Biomechanical Engineering, vol. 99, no. 3, pp. 126–147, 1977. [85] J. P. Garcia Jr, P. T. Garcia, and R. B. Rosen, “Retinal blood flow in the normal human eye using the canon laser blood flowmeter,” Ophthalmic research, vol. 34, no. 5, pp. 295–299, 2002. [86] C. E. Riva, J. E. Grunwald, S. H. Sinclair, and B. Petrig, “Blood velocity and volumetric flow rate in human retinal vessels.” Investigative ophthalmology & visual science, vol. 26, no. 8, pp. 1124–1132, 1985. [87] R. N. Weinreb and A. Harris, Ocular blood flow in glaucoma. Kugler Publications, 2009, vol. 6. 127 [88] V. Diliene, A. Harris, I. Januleviciene, L. Siaudvytyte, R. Barsauskaite, and B. Siesky, “Ocular hemodynamic changes after trabeculectomy in pseudoexfoliative and primary open-angle glaucoma,” Investigative Ophthalmology & Visual Science, vol. 54, no. 15, pp. 4445–4445, 2013. [89] S. Landahl, C. Bengtsson, J. A. Sigurdsson, A. Svanborg, and K. Svärdsudd, “Age-related changes in blood pressure.” Hypertension, vol. 8, no. 11, pp. 1044– 1049, 1986. [90] E. Pinto, “Blood pressure and ageing,” Postgraduate medical journal, vol. 83, no. 976, pp. 109–114, 2007. [91] S. S. Franklin, W. Gustin, N. D. Wong, M. G. Larson, M. A. Weber, W. B. Kannel, and D. Levy, “Hemodynamic patterns of age-related changes in blood pressure the framingham heart study,” Circulation, vol. 96, no. 1, pp. 308–315, 1997. [92] H. Qi, X. Bai, H. Zhou, and B. Wu, “Elasticity of blood vessel decreased induced by aging is the main factor of vascular senescence,” Heart, vol. 98, no. Suppl 2, pp. E146–E146, 2012. [93] D. A. Duprez, D. R. Jacobs Jr, P. L. Lutsey, D. Herrington, D. Prime, P. Ouyang, R. G. Barr, and D. A. Bluemke, “Race/ethnic and sex differences in large and small artery elasticity–results of the multi-ethnic study of atherosclerosis (mesa),” Ethnicity & disease, vol. 19, no. 3, p. 243, 2009. [94] D. W. Jones and J. E. Hall, “Racial and ethnic differences in blood pressure biology and sociology,” Circulation, vol. 114, no. 25, pp. 2757–2759, 2006. [95] B. Jani and C. Rajkumar, “Ageing and vascular ageing,” Postgraduate medical journal, vol. 82, no. 968, pp. 357–362, 2006. [96] K. J. Airaksinen, P. I. Salmela, M. K. Linnaluoto, M. J. Ikäheimo, K. Ahola, and L. J. Ryhänen, “Diminished arterial elasticity in diabetes: association with fluorescent advanced glycosylation end products in collagen,” Cardiovascular Research, vol. 27, no. 6, pp. 942–945, 1993. [97] M. A. Fazio, R. Grytz, J. S. Morris, L. Bruno, S. K. Gardiner, C. A. Girkin, and J. C. Downs, “Age-related changes in human peripapillary scleral strain,” Biomechanics and modeling in mechanobiology, vol. 13, no. 3, pp. 551–563, 2014. [98] M. J. Girard, J.-K. F. Suh, M. Bottlang, C. F. Burgoyne, and J. C. Downs, “Scleral biomechanics in the aging monkey eye,” Investigative ophthalmology & visual science, vol. 50, no. 11, pp. 5226–5237, 2009. [99] D. Yan, S. McPheeters, G. Johnson, U. Utzinger, and J. P. V. Geest, “Microstructural differences in the human posterior sclera as a function of age and race,” Investigative ophthalmology & visual science, vol. 52, no. 2, pp. 821–829, 2011. [100] R. Grytz, I. A. Sigal, J. W. Ruberti, G. Meschke, and J. C. Downs, “Lamina cribrosa thickening in early glaucoma predicted by a microstructure motivated growth and remodeling approach,” Mechanics of Materials, vol. 44, pp. 99–109, 2012. 128 [101] R. Grytz, C. A. Girkin, V. Libertiaux, and J. C. Downs, “Perspectives on biomechanical growth and remodeling mechanisms in glaucoma,” Mechanics research communications, vol. 42, pp. 92–106, 2012. [102] R. E. Norman, J. G. Flanagan, I. A. Sigal, S. M. Rausch, I. Tertinegg, and C. R. Ethier, “Finite element modeling of the human sclera: influence on optic nerve head biomechanics and connections with glaucoma,” Experimental eye research, vol. 93, no. 1, pp. 4–12, 2011. [103] I. A. Sigal, J. G. Flanagan, I. Tertinegg, and C. R. Ethier, “Finite element modeling of optic nerve head biomechanics,” Investigative ophthalmology & visual science, vol. 45, no. 12, pp. 4378–4387, 2004. [104] I. A. Sigal, J. G. Flanagan, and C. R. Ethier, “Factors influencing optic nerve head biomechanics,” Investigative ophthalmology & visual science, vol. 46, no. 11, pp. 4189–4199, 2005. [105] I. A. Sigal and C. R. Ethier, “Biomechanics of the optic nerve head,” Experimental eye research, vol. 88, no. 4, pp. 799–807, 2009. [106] I. A. Sigal, J. G. Flanagan, I. Tertinegg, and C. R. Ethier, “Modeling individualspecific human optic nerve head biomechanics. part i: Iop-induced deformations and influence of geometry,” Biomechanics and modeling in mechanobiology, vol. 8, no. 2, pp. 85–98, 2009. [107] ——, “Modeling individual-specific human optic nerve head biomechanics. part ii: influence of material properties,” Biomechanics and modeling in mechanobiology, vol. 8, no. 2, pp. 99–109, 2009. [108] I. A. Sigal, “An applet to estimate the iop-induced stress and strain within the optic nerve head,” Investigative ophthalmology & visual science, vol. 52, no. 8, pp. 5497–5506, 2011. [109] J. B. Jonas, “Ophthalmodynamometric assessment of the central retinal vein collapse pressure in eyes with retinal vein stasis or occlusion,” Graefe’s archive for clinical and experimental ophthalmology, vol. 241, no. 5, pp. 367–370, 2003. [110] J. Jonas, “Ophthalmodynamometric determination of the central retinal vessel collapse pressure correlated with systemic blood pressure,” British journal of ophthalmology, vol. 88, no. 4, pp. 501–504, 2004. [111] J. B. Jonas, “Ophthalmodynamometric measurement of orbital tissue pressure in thyroid-associated orbitopathy,” Acta Ophthalmologica Scandinavica, vol. 82, no. 2, pp. 239–239, 2004. [112] Y. Sugiura, F. Okamoto, Y. Okamoto, Y. Hasegawa, T. Hiraoka, and T. Oshika, “Ophthalmodynamometric pressure in eyes with proliferative diabetic retinopathy measured during pars plana vitrectomy,” American journal of ophthalmology, vol. 151, no. 4, pp. 624–629, 2011. [113] M. R. Glucksberg and R. Dunn, “Direct measurement of retinal microvascular pressures in the live, anesthetized cat,” Microvascular research, vol. 45, no. 2, pp. 158–165, 1993. 129 [114] R. Attariwala, C. P. Giebs, and M. R. Glucksberg, “The influence of elevated intraocular pressure on vascular pressures in the cat retina.” Investigative ophthalmology & visual science, vol. 35, no. 3, pp. 1019–1025, 1994. [115] S. J. Denardo, E. G. Yamada, V. K. Hargrave, and P. G. Yock, “Effect of stenosis inlet geometry on shear rates of blood flow in the upstream region,” American heart journal, vol. 125, no. 2, pp. 350–356, 1993. [116] C. R. Woodman, D. W. Trott, and M. H. Laughlin, “Short-term increases in intraluminal pressure reverse age-related decrements in endothelium-dependent dilation in soleus muscle feed arteries,” Journal of Applied Physiology, vol. 103, no. 4, pp. 1172–1179, 2007. [117] A. Boltz, D. Schmidl, R. M. Werkmeister, M. Lasta, S. Kaya, S. Palkovits, R. Told, K. J. Napora, A. Popa-Cherecheanu, G. Garhöfer et al., “Regulation of optic nerve head blood flow during combined changes in intraocular pressure and arterial blood pressure,” Journal of Cerebral Blood Flow & Metabolism, vol. 33, no. 12, pp. 1850–1856, 2013. [118] S. Cassani, A. Harris, B. Siesky, and J. Arciero, “Theoretical analysis of the relationship between changes in retinal blood flow and ocular perfusion pressure,” Journal of Coupled Systems and Multiscale Dynamics, vol. 3, no. 1, pp. 38–46, 2015. [119] S. Cassani, J. Arciero, G. Guidoboni, B. Siesky, and A. Harris, “Theoretical predictions of metabolic flow regulation in the retina,” Journal for Modeling in Ophthalmology, under review. [120] Z. He, J. K. Lim, C. T. Nguyen, A. J. Vingrys, and B. V. Bui, “Coupling blood flow and neural function in the retina: a model for homeostatic responses to ocular perfusion pressure challenge,” Physiological reports, vol. 1, no. 3, p. e00055, 2013. [121] T. Tani, T. Nagaoka, S. Nakabayashi, T. Yoshioka, and A. Yoshida, “Autoregulation of retinal blood flow in response to decreased ocular perfusion pressure in cats: Comparison of the effects of increased intraocular pressure and systemic hypotensionretinal blood flow in response to decreased opp,” Investigative ophthalmology & visual science, vol. 55, no. 1, pp. 360–367, 2014. [122] P.-y. Teng, J. Wanek, N. P. Blair, and M. Shahidi, “Response of inner retinal oxygen extraction fraction to light flicker under normoxia and hypoxia in ratinner retinal oxygen extraction fraction,” Investigative Ophthalmology & Visual Science, vol. 55, no. 9, pp. 6055–6058, 2014. [123] L. M. Haugh, R. A. Linsenmeier, and T. K. Goldstick, “Mathematical models of the spatial distribution of retinal oxygen tension and consumption, including changes upon illumination,” Annals of biomedical engineering, vol. 18, no. 1, pp. 19–36, 1990. [124] M. W. Roos, “Theoretical estimation of retinal oxygenation during retinal artery occlusion,” Physiological Measurement, vol. 25, no. 6, p. 1523, 2004. [125] S. J. Cringle and D.-Y. Yu, “A multi-layer model of retinal oxygen supply and consumption helps explain the muted rise in inner retinal po 2 during systemic hyperoxia,” Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology, vol. 132, no. 1, pp. 61–66, 2002. 130 [126] S. L. Polyak, “The retina: the anatomy and the histology of the retina in man, ape, and monkey, including the consideration of visual functions, the history of physiological optics, and the histological laboratory technique.” 1941. [127] P. Causin, G. Guidoboni, F. Malgaroli, R. Sacco, and A. Harris, “Blood flow mechanics and oxygen transport and delivery in the retinal microcirculation: multiscale mathematical modeling and numerical simulation,” Biomechanics and Modeling in Mechanobiology, pp. 1–18, 2015. [128] R. A. Linsenmeier and R. Braun, “Oxygen distribution and consumption in the cat retina during normoxia and hypoxemia.” The Journal of general physiology, vol. 99, no. 2, pp. 177–197, 1992. [129] B. W. Pogue, J. A. O’Hara, C. M. Wilmot, K. D. Paulsen, and H. M. Swartz, “Estimation of oxygen distribution in rif-1 tumors by diffusion model-based interpretation of pimonidazole hypoxia and eppendorf measurements,” Radiation Research, vol. 155, no. 1, pp. 15–25, 2001. [130] R. Ganfield, P. Nair, and W. Whalen, “Mass transfer, storage, and utilization of o2 in cat cerebral cortex,” American Journal of Physiology–Legacy Content, vol. 219, no. 3, pp. 814–821, 1970. [131] J. C. Arciero, B. E. Carlson, and T. W. Secomb, “Theoretical model of metabolic blood flow regulation: roles of atp release by red blood cells and conducted responses,” American Journal of Physiology-Heart and Circulatory Physiology, vol. 295, no. 4, pp. H1562–H1571, 2008. [132] D. Goldman and A. S. Popel, “A computational study of the effect of capillary network anastomoses and tortuosity on oxygen transport,” Journal of Theoretical Biology, vol. 206, no. 2, pp. 181–194, 2000. [133] Z. Li, T. Yipintsoi, and J. B. Bassingthwaighte, “Nonlinear model for capillarytissue oxygen transport and metabolism,” Annals of biomedical engineering, vol. 25, no. 4, pp. 604–619, 1997. [134] J. Moore and C. Ethier, “Oxygen mass transfer calculations in large arteries,” Journal of biomechanical engineering, vol. 119, no. 4, pp. 469–475, 1997. [135] P. Zunino, “Mathematical and numerical modeling of mass transfer in the vascular system,” Ph.D. dissertation, ÉCOLE POLYTECHNIQUE FÉDÉRALE DE LAUSANNE, 2002. [136] C. D’ANGELO, “Multiscale modelling of metabolism and transport phenomena in living tissues,” Ph.D. dissertation, ÉCOLE POLYTECHNIQUE FÉDÉRALE DE LAUSANNE, 2007. [137] R. A. Moses, “Intraocular pressure,” Adler’s Physiology of the Eye: clinical application. CV Mosby Co, pp. 223–245, 1987. [138] J. Kiel, “Physiology of the intraocular pressure,” in Pathophysiology of the Eye: Glaucoma, J. Feher, Ed. Budapest: Akademiai Kiado, 1998, pp. 109–144. [139] J. Kiel, M. Hollingsworth, R. Rao, M. Chen, and H. Reitsamer, “Ciliary blood flow and aqueous humor production,” Progress in retinal and eye research, vol. 30, no. 1, pp. 1–17, 2011. 131 [140] M. Szopos, S. Cassani, G. Guidoboni, C. Prud’homme, R. Sacco, B. Siesky, and A. Harris, “Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma treatment,” Journal for Modeling in Ophthalmology, in press. [141] G. Lyubimov, I. Moiseeva, and A. Stein, “Dynamics of the intraocular fluid: Mathematical model and its main consequences,” Fluid Dynamics, vol. 42, no. 5, pp. 684–694, 2007. [142] R. F. Brubaker, “The effect of intraocular pressure on conventional outflow resistance in the enucleated human eye.” Investigative ophthalmology & visual science, vol. 14, no. 4, pp. 286–292, 1975. [143] K. A. Johnson and R. S. Goody, “The original michaelis constant: translation of the 1913 michaelis–menten paper,” Biochemistry, vol. 50, no. 39, pp. 8264–8269, 2011. [144] I. Sobol, “Sensitivity analysis for nonlinear mathematical models,” Math Model Comput Exp, vol. 1, pp. 407–14, 1993. [145] G. Saporta, Probabilités, analyse des données et statistique. Editions Technip, 2006. [146] B. Sudret, “Global sensitivity analysis using polynomial chaos expansions,” Reliability Engineering & System Safety, vol. 93, no. 7, pp. 964–979, 2008. [147] R. S. Carel, A. D. Korczyn, M. Rock, and I. Goya, “Association between ocular pressure and certain health parameters,” Ophthalmology, vol. 91, no. 4, pp. 311–314, 1984. [148] H. D. Sesso, M. J. Stampfer, B. Rosner, C. H. Hennekens, J. M. Gaziano, J. E. Manson, and R. J. Glynn, “Systolic and diastolic blood pressure, pulse pressure, and mean arterial pressure as predictors of cardiovascular disease risk in men,” Hypertension, vol. 36, no. 5, pp. 801–807, 2000. [149] Y. H. Kwon, J. H. Fingert, M. H. Kuehn, and W. L. Alward, “Primary openangle glaucoma,” New England Journal of Medicine, vol. 360, no. 11, pp. 1113– 1124, 2009. [150] W. D. Stamer and T. S. Acott, “Current understanding of conventional outflow dysfunction in glaucoma,” Current opinion in ophthalmology, vol. 23, no. 2, p. 135, 2012. [151] A. H. Rulo, E. L. Greve, H. C. Geijssen, and P. F. Hoyng, “Reduction of intraocular pressure with treatment of latanoprost once daily in patients with normalpressure glaucoma,” Ophthalmology, vol. 103, no. 8, pp. 1276–1282, 1996. [152] J. H. Siggers and C. R. Ethier, “Fluid mechanics of the eye,” Annual Review of Fluid Mechanics, vol. 44, pp. 347–372, 2012. [153] E. A. Lippa, L.-E. Carlson, B. Ehinger, L.-O. Eriksson, K. Finnström, C. Holmin, S.-E. G. Nilsson, K. Nyman, C. Raitta, A. Ringvold et al., “Dose response and duration of action of dorzolamide, a topical carbonic anhydrase inhibitor,” Archives of ophthalmology, vol. 110, no. 4, pp. 495–499, 1992. 132 [154] A. Bayer, J. D. Henderer, T. Kwak, J. Myers, J. Fontanarosa, and G. L. Spaeth, “Clinical predictors of latanoprost treatment effect,” Journal of glaucoma, vol. 14, no. 4, pp. 260–263, 2005. [155] R. van der Valk, C. A. Webers, J. S. Schouten, M. P. Zeegers, F. Hendrikse, and M. H. Prins, “Intraocular pressure–lowering effects of all commonly used glaucoma drugs: a meta-analysis of randomized clinical trials,” Ophthalmology, vol. 112, no. 7, pp. 1177–1185, 2005. [156] P. Watson, J. Stjernschantz, L. S. Group et al., “A six-month, randomized, double-masked study comparing latanoprost with timolol in open-angle glaucoma and ocular hypertension,” Ophthalmology, vol. 103, no. 1, pp. 126–137, 1996. [157] N. Orzalesi, L. Rossetti, A. Bottoli, and P. Fogagnolo, “Comparison of the effects of latanoprost, travoprost, and bimatoprost on circadian intraocular pressure in patients with glaucoma or ocular hypertension,” Ophthalmology, vol. 113, no. 2, pp. 239–246, 2006. [158] R. Sacco, S. Cassani, G. Guidoboni, M. Szopos, C. Prud’homme, and A. Harris, “Modeling the coupled dynamics of ocular blood flow and production and drainage of aqueous humor,” Proceedings of the 4th International Conference on Computational and Mathematical Biomedical Engineering - CMBE2015, June 29-July 1, 2015, Cachan (France), pp. 608–611, 2015. VITA 133 VITA Simone Cassani Education • Doctor of Philosophy, August 2016 Purdue University, Applied Mathematics Blood circulation and aqueous humor flow in the eye: multi-scale modeling and clinical applications Major Advisor: Giovanna Guidoboni, IUPUI Co-advisors: Julia Arciero, IUPUI, Alon Harris, Indiana University • Master of Science, July 2011 Politecnico di Milano, Mathematical Engineering Implementation of Large Eddy Simulation models for turbulence on finite elements codes Advisors: Antonella Abbà and Lorenzo Valdettaro, Politecnico di Milano. • Bachelor of Science, March 2008 Politecnico di Milano, Mathematical Engineering Cooling and ice formation process in turbulent convection Advisors: Lorenzo Valdettaro and Antonella Abbà, Politecnico di Milano. Travel Awards • 2013, 2015, 2016 School of Science Graduate Student Council Travel Grant, IUPUI 134 • 2014, 2016 GPSG Educational Enhancement Travel Grant, IUPUI • 2013 Graduate Student Travel Fellowship Award, IUPUI Teaching Experiences • Recitation instructor MATH 16600: Integrated Calculus and Analytic Geometry II, Fall 2014, Primary instructor J. Miller, IUPUI. • Primary instructor MATH 15400: Algebra and Trigonometry II, Spring 2014 and Fall 2014, Coordinator S. Rangazas, IUPUI. • Primary instructor MATH 15300: Algebra and Trigonometry I, Fall 2013, Coordinator S. Rangazas, IUPUI. • Recitation instructor MATH 16600: Integrated Calculus and Analytic Geometry II, Spring 2013, Primary instructor J. Watt, IUPUI. • Tutor at the Mathematics Assistance Center (MAC), IUPUI, 2013-2016. Research Experiences • January-December 2015, Université de Strasbourg, Labex IRMA, France. Research project: Eye2Brain: the eye as a window on the brain. Supervisors: Giovanna Guidoboni, IUPUI, Christophe Prud’homme, IRMA and Marcela Szopos, IRMA. • March 2014 and 2016, SQuaREs (Structured Quartet Research Ensembles), American Institute of Mathematics. Research project: Ocular blood flow and its role in development of glaucoma. Sergey Lapin, Giovanna Guidoboni, Julia Arciero, Alon Harris, Lucia Carichino, Simone Cassani. • August-December 2011, Short-term Scholar, IUPUI. Research project: Mathematical modeling of the ocular blood flow. Supervisors: Giovanna Guidoboni, IUPUI, Alon Harris, IU School of Medicine. 135 Conferences • The ARVO Annual Meetings, May 2012, 2013, 2014, 2015 and 2016, USA. • The International Congress of Advanced Technologies and Treatments for Glaucoma, October 2015, Milano, Italy. • 11th World Congress on Computational Mechanics, July 2014, Barcelona, Spain. • Integrated Multidisciplinary approaches in the study and Care of the Human Eye, Clinical Experience - Mathematical Modeling - New Technologies, June 2013, Milano, Italy. • 13th INFORMS Computing Society Conference, January 2013, Santa Fe, NM, USA. • The Third and Fourth Annual Eugene & Marilyn Glick Eye Institute Vision Symposium, October 2013 and November 2012, Indianapolis, IN, USA. Selected Publications • Theoretical predictions of metabolic flow regulation in the retina, S. Cassani, J. Arciero, G. Guidoboni, B. Siesky, A. Harris, Journal for Modeling in Ophthalmology, Accepted (May 2016) • Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management, M. Szopos, S. Cassani, G. Guidoboni, C. Prud’homme, R. Sacco, B. Siesky, A. Harris, Journal for Modeling in Ophthalmology, Accepted (April 2016) • Intraocular Pressure, blood pressure and retinal blood flow autoregulation: a mathematical model to clarify their relationship and clinical relevance, G. Guidoboni, A. Harris, S. Cassani, J. Arciero, B. Siesky, A. Amireskandari, L. Tobe, P. Egan, I. Januleviciene, J. Park, Invest Ophthalmol Vis Sci, 55(7) (2014) 4105-18 136 • Effect of trabeculectomy on retinal hemodynamics: mathematical modeling of clinical data, S. Cassani, G. Guidoboni, I. Januleviciene, L. Carichino, B.A. Siesky, L.A. Tobe, A. Amireskandari, D.P. Baikstiene, A. Harris, Integrated Multidisciplinary Approaches in the Study and Care of the Human Eye, Kugler Publications, (2014) 29-36 • Mathematical modeling approaches in the study of glaucoma disparities among people of African and European descents, G. Guidoboni, A. Harris, J.C. Arciero, B.A. Siesky, A. Amireskandari, A.L. Gerber, A.H. Huck, N.J. Kim, S. Cassani and L. Carichino, Journal of Coupled Systems and Multiscale Dynamics, Vol. 1(1) (2013) 1-21 Other Academic Activities • 2012-2016 Officer of the Society for Industrial and Applied Mathematics (SIAM) Student Chapter at IUPUI • 2013 International Graduate Welcome Volunteer, IUPUI