UNIVERSITÀ DEGLI STUDI DI MILANO

Facoltà di Scienze Matematiche, Fisiche e Naturali

Corso di Laurea Magistrale in Matematica

MULTISCALEMODELS AND

NUMERICALSIMULATIONOFRETINALMICROCIRCULATION:

BLOODFLOW AND

MASSTRANSPORTPHENOMENA

Relatore: Dr.ssa Paola CAUSINCorrelatore: Prof. Riccardo SACCO

Tesi di Laurea di:Francesca MALGAROLIMatricola n. 791674

Anno Accademico 2011/2012

Contents

1 Models of Retinal Geometry 6

1.1 Anatomical overview of the retina. . . . . . . . . . . . . . . . 61.1.1 Retinal vascular network. . . . . . . . . . . . . . . . . 8

1.2 Geometrical Models of the Retina . . . . . . . . . . . . . . . . 131.2.1 Arteriolar Tree . . . . . . . . . . . . . . . . . . . . . . 13

1.2.1.1 Dichotomic Tree . . . . . . . . . . . . . . . . 131.2.1.2 Diusion-limited aggregation model (DLA) . 16

1.2.2 Capillary Plexi . . . . . . . . . . . . . . . . . . . . . . 191.2.3 Tissue with layers . . . . . . . . . . . . . . . . . . . . 21

2 Models of Blood Flow 23

2.1 Flow in the arteriolar tree . . . . . . . . . . . . . . . . . . . . 232.1.1 Model equations . . . . . . . . . . . . . . . . . . . . . 262.1.2 Numerical results . . . . . . . . . . . . . . . . . . . . . 30

2.2 Flow in Capillary Plexi and in Interstitial Tissue . . . . . . . 442.2.1 Calculation of Permeability in Capillary Bed . . . . . 45

3 Models for Solute Transport and Delivery 47

3.1 Solute transport across the vessel wall . . . . . . . . . . . . . 473.1.1 The Wall-Free Model . . . . . . . . . . . . . . . . . . . 493.1.2 The Multilayer Model . . . . . . . . . . . . . . . . . . 513.1.3 The case of solute oxygen: free oxygen and oxygen

carried by oxyhemoglobin . . . . . . . . . . . . . . . . 543.2 Solute Transport in Capillary Beds and in Interstitial Tissue . 573.3 Solute transport in the tissue . . . . . . . . . . . . . . . . . . 58

3.3.1 The O2 model for Retinal Tissue . . . . . . . . . . . . 603.3.2 Numerical Solution . . . . . . . . . . . . . . . . . . . . 60

4 Multiscale coupled Model 64

4.1 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . 66

A Calculation of Equivalent Resistance in Capillaries Bed 77

B The Finite Element Method 80

C The Scharfetter-Gummel method 85

Bibliography 88

1

Riassunto

La retina è il solo tessuto dell'organismo nel quale i vasi sanguigni possonoessere studiati in vivo in maniera non invasiva. Tuttavia, i numerosi processiinerenti la microcircolazione che vi hanno luogo sono complessi e non an-cora completamente compresi. La circolazione oculare svolge infatti funzionidelicate, essendo capace di reagire a numerosi stimoli dierenti e mantenereuna condizione di omeostasi. Lo studio di tali meccanismi regolatori in con-dizioni siologiche è fondamentale per poter porre in atto terapie adeguatenel caso in cui essi vengano a mancare, causando gravi patologie quali adesempio il glaucoma, seconda causa di cecità nel mondo. A questo scopo, imodelli matematici e computazionali basati sui principi della meccanica, ui-dodinamica, trasporto di massa e elettrochimica, possono fornire indicazionipreziose per comprendere i processi che hanno luogo nella microcircolazioneretinale. In questo lavoro di tesi, viene proposto un modello matematico perlo studio della microcircolazione retinale nei suoi dierenti compartimenti(arteriole, letti capillari e tessuto). Viene dapprima costruita una strutturageometrica articiale attraverso un algoritmo di Diusion Limited Aggrega-tion (DLA) che rappresenti in modo adeguato la complessa architettura deivasi sanguigni maggiori. Successivamente, vengono studiati i modelli per lacircolazione sanguigna e per la diusione di soluti, con particolare riferimentoalla dinamica dell'ossigeno. Tali modelli rappresentano i distretti in esame alivello microscopico (tessuto), mesoscopico (arteriole) o macroscopico (letticapillari). Vengono condotte simulazioni numeriche basate su risolutori adelementi niti, ottenendo i proli di velocità del usso sanguigno nei varidistretti, che vengono validati sulla base di dati sperimentali di letteratura.Viene inoltre calcolata negli stessi distretti la pressione parziale di ossigeno,il cui valore - se non mantenuto entro valori siologici - può condurre alleretinopatie di cui accennato sopra.

La tesi è organizzata come segue:

• nel Capitolo 1, viene fornita una breve descrizione anatomica dellaretina con i suoi distretti e si discute la relativa rappresentazione geo-metrica usata nel modello matematico (albero dicotomico oppure DLA

2

per la rete di arteriole, geometria regolarizzata per i letti capillari,struttura a strati per il tessuto retinale);

• nel Capitolo 2, vengono presentati i modelli matematici per lo studiodel usso sanguigno nella rete di arteriole e nei letti capillari. Nel primocaso, viene adottato un prolo di velocità di tipo Hagen-Poiseuille, stu-diato nella reticolazione tramite leggi di conservazione. I letti capillarivengono invece rappresentati con un modello di tipo 0D, tramite unaresistenza equivalente calcolata sulla base della loro geometria sem-plicata. Tale resistenza viene accoppiata alla rete arteriolare, costi-tuendone il carico nale. Viene inoltre discusso un possibile modellopiù dettagliato per la rappresentazione del letto capillare tramite unapproccio di tipo Double Continuum, che descrive il usso nei vasi ca-pillari e il usso interstiziale nel tessuto circostante come due processidivisi ma accoppiati tramite leggi di scambio di massa;

• nel Capitolo 3, viene studiata la dinamica del trasporto di soluti nei di-stretti retinali, incluso il tessuto in cui risiedono i neuroni fotorecettori,con particolare riferimento all'ossigeno. In particolare, viene studiatoil trasporto di ossigeno, tramite modelli di tipo Wall-Free oppure Mul-tilayer, anche attraverso la parete dei vasi della rete arteriolare, doveil suo livello negli strati di endotelio e di cellule muscolari lisce è notoessere responsabile - tramite una complessa catena di fenomeni chimici- della regolazione del diametro dei vasi stessi;

• nel Capitolo 4, vengono discusse le metodologie di accoppiamento nu-merico dei vari modelli coinvolti (usso sanguigno e trasporto di soluto)nei vari distretti (rete arteriolare, letti capillari, tessuto). Vengono inol-tre proposti alcuni temi che si ritengono interessanti da sviluppare inun futuro lavoro di ricerca su questo argomento.

3

Introduction

Retina is the only tissue in which blood vessels can be studied in vivo in anon-invasive way, the many processes that take place and their relationshipsbeing numerous and dicult to study. Ocular circulation is a delicate mech-anism, charged to maintain the homeostasis of retinal function in responseto physiological stimuli. It is thus crucial to understand the processes under-lying the regulation of ocular circulation in physiological conditions. Theirimpairment causes severe retinal disorders, like glaucoma that is the secondcause of blindness worldwide [13].Mathematical and computational models based on the physical principlesof mechanics, uid-dynamics, mass transport and electrochemistry can helpunraveling the cause-eect mechanisms acting as key factors in the regulationof functioning of the retinal microvasculature [8].

In this thesis, we propose mathematical models for the study of the abovementioned phenomena in the dierent districts of the retina (arterioles, cap-illaries and tissue). We also propose an articial geometric structure torepresent the complex network of blood vessels in the plane of the arterioles.We have divided the work into three parts: in the rst one we present thegeometric models used to describe the retina, in the second the mathemat-ical models for the blood ow and in the last the mathematical models formass transport.

In the rst chapter we give a brief anatomical description of the retinaand we present the geometrical model used to describe it. To describe thearterial network, we present two models. The rst is an existing model char-acterized by a dichotomic symmetric branching system. The second is amodel that is based on the Diusion-Limited Aggregation model (DLA). Wechose this model because numerous studies have shown that the fractal di-mension of the structures that are obtained with this approach is very similarto the fractal dimension observed from images of the retinal fundus.The retinal capillary plexus is formed by capillaries that are embedded inthe retinal tissue and each of these have to be described by mathemati-cal equations very dierent from each other. Therefore we must treat thetissue and the capillaries as two separated domains. To solve this prob-

4

lem we use a geometric model derived from petroleum reservoir analysis,the double-continuum approach, in which the two continua are formed bytissue-interstitial uid and tissue-capillaries.Finally, the retinal tissue, for its structure, is treated as a multilayered do-main divided into eight levels each of which has dierent characteristics.

In the second chapter, we present the mathematical models for bloodow in the arterial tree and in the capillary plexus. In the arterial tree,blood ow is considered as a Poiseuille ow and we highlight the similaritybetween geometric models presented and an electrical circuit. Then, assum-ing in each node of bifurcation the mass conservation law for uid ow, theresulting system of equations gives the mean hemodynamic parameters ofthe network in each vessel. To treat the coupling between the arterial treeand the capillary plane, we insert in the system an equivalent resistance thatwe calculate considering the capillary plane as a circuit formed from squaremeshes. This is presented in appendix. Finally, the numerical solutions arepresented and are compared with the experimental data.

In the third chapter we study the solute transport within the retina. Atthe beginning we present general mathematical models which are then spe-cialized for oxygen. We present two models that study the solute transportin arterial tree considering that the solute, in addition to being transportedand diused along the vessel, lters also out through the wall of the vesselitself. The vessel wall can be considered by a "transfer" boundary condition(Wall-Free model) or formed by three layers (endothelium, smooth musclecells layer and tissue) with dierent characteristics from each other (Multi-layer model).Within the tissue, we nd the plane of the capillaries that provides oxygento the tissue itself and the planes in which are found the largest number ofphotoreceptors, and therefore the greater consumption of oxygen, which com-municates with the brain to generate the visual image. For this, the oxygentransport in the tissue is described by a system of diusion equations withsource and consumption terms. These terms depend on the characteristicsof the layer that is considered.

In the fourth and nal chapter, we discuss the coupling between all mod-els we presented in the previous chapters, both for blood ow and solutetransport, in each district of the retina (arterial tree, capillary beds and tis-sue). The numerical results are discussed. We also propose some possiblethemes that can be developed for future research work on this topic.

5

Chapter 1

Models of Retinal Geometry

Awareness of the uniqueness of the retinal vascular patterns dates back to1935 when two ophthalmologists, C. Simon and I.Goldstein [18], while study-ing eye diseases, realized that every eye has its own totally unique patternof blood vessels. P.Tower, studying identical twins [21], showed in his studythat of all the factors compared between the twins, retinal vascular patternsshowed the least similarities. We are thus faced with a problem in a complexgeometry, strongly varying among individuals. In the following, we providea brief description of retinal anatomy and we discuss the geometrical modelswe adopt for its mathematical study.

1.1 Anatomical overview of the retina.

In adult humans, the retina is approximately 72% of a sphere about 22mm indiameter and lines the back of the eye. The retina is a light-sensitive tissuelining the inner surface of the eye. The optics of the eye create an image ofthe visual world on the retina, which serves much the same function as thelm in a camera. Light striking the retina initiates a cascade of chemicaland electrical events that ultimately trigger nerve impulses. These are sentto various visual centres of the brain through the bres of the optic nerve. Invertebrate embryonic development, the retina and the optic nerve originateas outgrowths of the developing brain, so the retina is considered part of thecentral nervous system (CNS) and is actually brain tissue. That is why oftentheir vasculature and hemodynamic parameters are treated as if they werethe same. A section of a portion of the retina reveals that is composed ofseveral layers (see Fig.1.1) in which ganglion cells (the output neurons of theretina) lie innermost closest to the lens and front of the eye while photore-ceptors (rods and cones) lie outermost against the pigment epithelium andchoroid. Light passes through several transparent nerve layers to reach therods and cones. A chemical change in the rods and cones send an electricalsignal back to the nerves. The signal goes rst to the bipolar and horizontal

6

cells, then to the amacrine cells and ganglion cells, then to the optic nervebres.

Figure 1.1: Section of the retina through its thickness. The image highlightsthe neural structures.

The optic disc, a part of the retina sometimes called "the blind spot"because it lacks photoreceptors, is located in the optic papilla, a nasal zonewhere the optic-nerve bres leave the eye. The optic disc appears as anoval white area of about 2x1.5mm. Temporal to this disc is the macula (seeFig. 1.2). At its center, approximately 15mm to the left of the disc, is thefovea, a spot measuring less than a quarter of a millimeter (200 µm) that isresponsible for sharp central vision. The circular eld of 6mm around thefovea is considered the central retina, while the area beyond this is calledthe peripheral retina.

The retinal thickness shows great variations (see Fig. 1.3). The retina isthinnest at the foveal oor (0.10, 0.150-0.200 mm) and thickest (0.23, 0.320mm) at the foveal rim. Beyond the fovea the retina rapidly thins until theequator. At the ora serrata the retina is thinnest (0.080 mm).

7

Figure 1.2: Sectional detail of the retina along the superior-inferior axis ofa left human eye through the optic nerve, showing details of the vascularsupply in this location.

Figure 1.3: The variation of human retinal thickness inmm around the fovea(data from Sigelman and Ozanics (1982)).

1.1.1 Retinal vascular network.

There are two sources of blood supply to the mammalian retina: the centralretina artery (CRA) and the choroidal blood vessels. The choroid vesselsprovide the greatest blood ow (65-85%) and are vital for the maintenanceof the outer retina (particularly the photoreceptores). The remaining 20-30%blood ow to the retina comes from the central retina artery (see Fig. 1.4).Within the optic nerve, the CRA divides to form two major trunks and eachof these divides again to form the superior nasal and temporal and the inferiornasal and temporal arteries that supply the four quadrants of the retina.Theretinal venous branches are distributed in a similar fashion. The vesselsemerge from the optic nerve and run in a radial fashion curving toward andaround the fovea. The major arterial and venous branches and the successive

8

divisions of the retinal vasculature are present in the nerve ber layer close tothe internal limiting membrane. The retinal arterial circulation in the humaneye is a terminal system with no arteriovenous anastomoses (communicationbetween vessels) or communication with other arterial systems: thus, theblood supply to a specic retinal quadrant comes exclusively from the specicretinal arteries and veins that supply that quadrant, and any blockage ofblood supply results into infarction.

Figure 1.4: Section of the retina with the distribution of blood vessels. Figuretaken from [2]

As the large arteries extend in the retina towards the periphery theydivide to form arteries with progressively smaller diameter, until they reachthe point where they return continuously to the venous drainage system.This process of division occurs either dichotomously or at right angles tothe original vessels. The terminal arterioles and venules form an extensivecapillary network in the inner retina as far as the external border of the innernuclear layer.The retinal vasculature is structured in three distinct layers: the supercial(innermost) layer, the intermediate layer, within ganglion cells layer, andthe deep layer, within the inner nuclear layer (see Fig. 1.4). The largervessels lay in the innermost layer, whereas a plexus of capillaries occupythe other two layers with precapillary arterioles and postcapillary venuleslinking them to larger vessels. Blood ow of the arterioles in the supercial

9

layer is directed to the intermediate and deep layers of the retina. Normally,no blood vessels from the the CRA extend into the outer plexiform layer,the layer that divides the photoreceptores layer to the other layers. Thus,the photoreceptor layer of the retina is free of the blood vessels supplied bythe CRA. The choriocapillaries provide the blood supply to photoreceptors.Since the fovea contains only photoreceptores, this cone-rich area is freeof any branches from the CRA. The walls of all blood vessels except thecapillaries are composed of three distinct layers, or tunics (see Fig. 1.5).The tunics surround a central blood-containing space called the lumen. Theinner most tunic, which is in intimate contact with the blood, is the tunicaintima. It contains the endothelium that lines the lumen of the vessel, and itsat cells t closely together, forming a slick surface that minimizes friction.The endothelial cells of the retinal arteries are linked by tight junctions thatestablish a blood-retinal barrier to prevent the movement of large moleculesor plasma proteins in or out of the retinal vessels. The middle layer, thetunica media consists mostly of one or more layer of smooth muscle cells(SMCs). Since small changes in blood vessel diameter greatly inuence bloodow and blood pressure, the activities of the tunica media are critical inregulating circulatory dynamics. The tunica media is usually the bulkiestlayer in arteries.

Figure 1.5: Schematic structure of the wall of blood vessels.

The outer most layer is the tunica adventitia and is composed of looselywoven collagen brils that protect the vessel and anchor it to surroundingstructures. Capillaries are the smallest blood vessels. Their exceedingly thinwalls consist of just a thin tunica intima. Unlike the arteries and veins, capil-laries are very thin and fragile. They are so thin that blood cells can only pass

10

through them in single le. The capillary wall is composed of three distinctelements: endothelial cells, intramural pericytes and a basement lamina. Theexchange of oxygen and carbon dioxide takes place through the thin capil-lary wall. The red blood cells inside the capillary release their oxygen whichpasses through the wall and into the surrounding tissue. The tissue releasesits waste products, like carbon dioxide, which passes through the wall andinto the red blood cells. The continuous endothelial cell layer is surroundedby a basal lamina within which pericytes form a discontinuous layer in al-most a one to one ratio with endothelial cells. Like SMCs, pericytes providestructural support to the vasculature and represent the myogenic mechanismfor vasculature autoregulation of blood ow in response to changes in neuralactivity. They are able to regulate the capillary diameter through contrac-tion and relaxation. Pericytes are also involved in the regulation of vascularpermeability.

Arteries and veins physiology. In the human retina arteries and veinsaccompany each other, but they are distinguished based on the branchingpattern and the size of the vessels. Pattern. The arteries tend to have'Y-shape' branches with arms of equal diameter at the equator and at theperiphery of the retina. They give rise to side-arm branches which then pro-gressively divide into dichotomic branches of arterioles. The arterioles giverise to capillaries. As in the arteries, side-arm branches also arise from theveins, and give rise to venule branches. Unlike the arterioles, the venules aremore likely to have a 'conveying type' branching pattern, which is also knownas strictly asymmetric branching. There are veins which are sensitively big-ger than other, and they have a 'T-shape' and give uneven size to the wholeveins network. The arterioles with the delivery branching pattern are morespaced out in comparison with the vessels with the conveting branching pat-tern. Measurements. The arteries around the optic nerve are approximately100µm in diameter, with 18µm thick walls - then they decrease in diameter,until the branched arteries lying in the deeper retina reach 15µm. The majorbranches of the central veins close to the optic disk have a lumen of nearly200µm with a thin wall made up of a single layer of endothelial cells havinga thin basement membrane (0.1µm). The lack of smooth muscle cells in thevenular vessel wall results in a loss of a rigid structural framework for suchvessels, resulting in shape changes under condition of sluggish blood ow(e.g, diabetes) or with increased venous pressure. The retinal arteries havea thicker muscular layer, which allows increased constriction in response topressure and chemical stimuli.Blood ow of the supercial layer containing large vessels is mostly directedto the vasculature at the intermediate and deep layers of the retina, as it isshown in Figure 1.6.

11

Capillary physiology. Pattern. The retinal capillary network is spreadthroughout the retina, diusely distributed between the arterial and venoussystems in the intermediate layer and in the deep layer, and it is anosto-motic. The capillaries are connected in tri-junction connection pattern, inwhich each capillary is connected to two other capillaries. The capillarieseither form a 'loop' shape if they are distributed in the same layer, or movetransversely to connect vessels in the other layers. In the human retina aregional variation of the density of capillary distribution is reported: the cap-illary distribution at the equator region is denser than that in the peripheralregion. Measurements. There are three specic areas of the retina that aredevoid of capillaries, the 400 µm wide area centred around the fovea, theone adjacent to the major vessels and the retinal periphery. The capillariesnetwork extends as far peripherally as retinal arteries and veins. The retinalcapillary lumen is extremely small (3.5-6 µm in diameter). Like capillarynetworks elsewhere in the body, the retinal capillaries assume a meshworkconguration to ensure adequate perfusion to all retinal cells. The deep cap-illary layer has mesh diameter (i.e., the distance between capillaries) thataverages 50µm in diameter but varies between 15 and 130 µm. The moresupercial capillary layer has slightly larger meshwork, on average 65µm indiameter (16 to 150µm). In the mid-equatorial and anterior zones, wherethe retina is thinner, only one capillary layer is present.

Figure 1.6: Schematic representation of the retinal layers and the connectionbetween vessels.

12

1.2 Geometrical Models of the Retina

Although it would be necessary to study uid ow and transport of solutes ina real structure of the retina from medical imaging, this is beyond the scopesof the present work. Moreover, we would need a large number of parameters,which in some cases are dicult to determine experimentally (for examplethe distribution of the large number of capillaries in the capillary beds).Therefore, we use articial mathematical structures to describe the retina.The results we get are in good agrement with experimental data showingthat we use an acceptable modelization of the real geometry. In this chapterwe discuss the geometrical models used to describe the arteriolar tree, thecapillary bed and the tissue, respectively.

1.2.1 Arteriolar Tree

We present two types of vascular tree structures that we use in our studies.First, we describe the structure proposed in [20], a dichotomic branchingsystem, then we describe a structure more similar to the real retina basedon fractals.

1.2.1.1 Dichotomic Tree

In their paper [20], Takahashi and Nagaoka develop a theoretical and math-ematical concept to quantitatively describe the hemodynamic behaviour inthe microvascular network of the human retina. A dichotomic symmetricbranching network of the retinal vasculature is constructed, based on a com-bination of Murray's law and a mathematical model of fractal vascular trees.The optimal branching structure of a vascular tree is given by Q = krm,where Q is the volumetric ow rate, r is the inner radius of the vessel seg-ment, k is a constant, and m a junction exponent which ranges between2.7 and 3, as shown [19][17]. In [20], the constant m is set equal to 2.85,a value that is more suitable for application to the retinal microcirculation(see Fig. 1.7). It is proved mathematically that the exponent m is the sum ofa fractal dimension (D) and a branch exponent (α). Takahashi et al. applya fractal dimension of 1.70 and a branch exponent of 1.15.The fractal dimension can quantify the property of a complex vascular net-work, and the value of α can also quantitatively dene the relation betweenthe length and radius of a branch segment as

L(r) = 7.4rα. (1.1)

The equation is derived from data on cerebral vessels, but it is known fromstudies that the vasculature of the retina and brain are similar.

13

Figure 1.7: Ratio of larger daughter-branch diameters to their mother-branchdiameters vs asymmetry ratio of the larger to the smaller daughter branchdiameters at some bifurcations in the human retina. Dotted and solid line:curves predicted by Murray's Law with diameter exponent 3 and 2.85 re-spectively. Scattered data from photographed normal human eye.(Figure takenfrom [20])

Combining Q = kr2.85 with conservation of mass ow rate, the congura-tion of a dichotomic vascular tree at every branching point can be expressedby

r2.851 = r2.85

2,1 + r2.852,2 (1.2)

where r1 is the radius of a mother branch, and r2,1 and r2,2 are the radii ofdaughter branches at the same bifurcation. The larger arteriole, that orig-inates directly from the central retinal artery (CRA), is given a generationnumber of 1. Branches of respectively arteriole 1 are given a generation num-ber of 2, and subsequent generations are formed in an identical fashion untilthe ospring decreases about 6 µm in diameter. Individual precapillary ves-sels instead spread out into four true capillaries vessels, and then join againto form a single postcapillary venule, as shown in Fig. 1.8.

14

Figure 1.8: Microvascular arterial network topologically represented as a suc-cessively repeating dichotomic branching system. Each parent vessel givesrise to two osprings, each of the osprings gives rise to further two o-springs, and so on. Four capillaries are assumed to divide from each precap-illary. (Figure taken from [20])

In our study, we consider a network with 10 generations. The diametersdecrease from the larger arteriole, which has a diameter of 108 µm, valuetaken as an input data for the model, through small arterioles to precapil-laries, with diameter of 12.1µm.

Figure 1.9: Distribution of diameters as a function of the hierarchical level.

15

1.2.1.2 Diusion-limited aggregation model (DLA)

The application of fractals and fractal growth processes to the branchingblood vessels of the normal human retinal circulation was introduced byMasters and Patt in 1989 [12]. Growing, branching objects can be repro-duced by computer simulations in which the spatial dependence of a eldsatises the Laplace equation with moving boundary conditions. A class ofprocesses based on fractal growth is the diusion-limited aggregation model(DLA) proposed by Witten and Sander in 1981 [23], and applicable to ag-gregation in any system where diusion is the primary means of transportin the system. DLA can be described as the process whereby a single par-ticle performs a random walk until it accidently hits an existing immobileaggregate. Then, the particle attaches to the cluster and becomes immobile.Besides this irreversible sticking no interaction of particles is present. Thisextremely simple process produces surprisingly complex, branched objectswhich are very appealing from a scientic point of view and have evidencein natural phenomena. This is, for instance, approximately the case whenwater molecules form a snow ake. Under appropriate conditions, the vicin-ity of the ake contains almost no water and molecules have to cross thiswater poor region by means of Brownian motion before they can attach tothe growing aggregate. Another example is the electrolytic deposition ofmaterial on an electrode. If the present electric eld is not too large, themotion of ions will be dominated by diusion.The simulation of a DLA yields branching patterns similar to the branchingpatterns seen in the human retina. Moreover, a series of papers led to anestimate of the fractal dimension for the retinal vessels of D = 1.7, whichis in good agreement with the dimension of a diusion-limited aggregationcluster grown in two dimensions that is usually D = 1.71 . The fractal di-mension is dened as D = ln(M)/ln(R) where R is the radius of the cluster,the maximum distance between the rst particle to another, andM its mass(i.e. the number of the particles that form the cluster).

In our study, we have used a Matlab code to obtain DLA clusters, anexample of which is depicted in Fig. 1.10. To build it, we consider a massM = 3000 and a maximum radius R = 300. Then we use the Hit-and-Misslter, a binary morphological operators, with few structuring elements tomodify the cluster (see Fig. 1.11), so that:

• segments of the cluster are formed by a single pixel in width (a proce-dure similar to skeletonization);

• segments do not intersect forming closed regions (loops);

• each line has at most two forks.

16

We note that though the mass of B in Fig. 1.10 is lower than the originalDLA, the fractal dimension is still similar to the value of retinal fractaldimension.Since the retina has an area of about 1089 mm2 (72% of a sphere about22mm in diameter), we set a pixel to be equivalent to approximately κ µm(κ depends on the size of the matrix representing the gure, in our caseκ= 162µm). In this way, we can determine the length of the various segmentsin the DLA system.

Figure 1.10: A) The original DLA cluster. B) The DLA cluster after theapplication of the morphological operator.

Figure 1.11: A) A particular of the original cluster. B), C) and D) are thesame object as in A) after some application of the morphological operator.In B) each segment has a width of one pixel, C) shows a deleted loop, whileD) shows a deleted trifurcation.

17

Then we determine the four major arteries, and we assign to these aradius that decreases as one moves away from the center, between 54 and15 µm. We assign to smaller branches a value which decreases moving awayfrom the main branches from 15 to 5 µm.

Figure 1.12: Colours represent the values of the radii of the vessel branches.Larger radii correspond to the four major arteries.

18

1.2.2 Capillary Plexi

The capillaries are the main location where the transport of nutrients be-tween blood and tissue takes place.The number of capillaries is huge and the structure of the capillary plexusis very complex. The retinal capillaries are embedded in the retinal tissueand each of these has to be described by mathematical laws very dierentfrom each other. Therefore we must treat the tissue and the capillaries astwo separated domains.In the retinal tissue its individual components are not densely packed. There-fore the interstitial uid can ow freely within the tissue. For this reason,the retinal tissue can be described as a porous medium with a single phaseow, where with phase we mean a matter that has a homogeneous chemicalcomposition and physical state (in this case uid). The composition of in-terstitial uid is similar to blood plasma, which consists by 90% of water.In the capillary bed, to avoid the high computational expenses that we in-cur if we use a discrete approach, can be introduced a capillary continuum,which represents the capillary bed around the tissue as averaged quantity.Hence, capillary bed can be described as a porous medium with a uid phase,the blood, where the medium is represented by the tissue. To pass from adiscrete to a continuum description, we use the concept of a representativeelementary volume (REV) and we dene new eective parameters, as poros-ity, tortuosity or permeability.For these reasons, we use a double continuum approach [5], in which the het-erogeneous domain of capillary plexus is represented by two separate, butspatially overlapping and interacting continua, one consisting of the capil-laries and the tissue and the other of the interstitial uid and the tissue(see Fig.1.13). At any point of the capillary plexus domain two values foreach eective parameter are dened: one for the capillaries and one for theinterstitial uid within tissue.

19

Figure 1.13: Schematic illustration of the general concept of the double-continuum model.(Figure taken from [5])

20

1.2.3 Tissue with layers

Anatomically, the retina is usually considered to consist of two main parts,the outer retina (which is avascular) and the inner retina (which is sup-plied with blood). Moreover, it has a distinctly layered structure in whichoxygen sources and consumption are more compartmentalised than in othertissues. The oxygen required in the retina is primarily derived from bloodin choroidal vessels and in the central retina artery. The choroidal bloodvessels supply oxygen to the outer retina whereas the central retina arterysupplies the inner retina.In this study, we assume that the retinal tissue is divided in eight layers(see Fig. 1.14), based on their anatomical and functional properties [3]-[4].The outer retina is formed by three layers: outer segments of the photore-ceptors layer, inner segments of the photoreceptors layer and outer nuclearlayer. The inner retina is divided in ve layers: outer plexiform layer, in-ner nuclear layer, outer region of the inner plexiform layer, inner region ofthe inner plexiform layer and Ganglion cell/nerve bre layer. Most of theoxygen delivered by choroidal circulation to the outer retina is consumed byinner segments of the photoreceptors because in this layer are localized themajority of photoreceptors. A greater portion of the oxygen provided bythe retinal circulation to the inner retina is utilized by the inner plexiformlayer. Therefore, we will assume that the consumption of oxygen takes placeonly in these two layers. In the Ganglion cell/nerve bre layer is locatedthe supercial capillary bed while in the outer plexiform layer is it the deepcapillary bed.

21

Figure 1.14: Scheme of retinal layers. Layer 1=outer segments of thepho-toreceptors layer, Layer 2=inner segments of the photoreceptors layer,Layer3 = and outer nuclear layer, Layer 4 = outer plexiform layer, Layer 5 = innernuclear layer, Layer 6 = outer region of the inner plexiform layer,Layer 7 =inner region of the inner plexiform layer, Layer 8 = Ganglion cell/nerve brelayer.(Figure taken from [3])

22

Chapter 2

Models of Blood Flow

The circulatory system in general, and so also the human retina we areanalyzing, consists of a network of many interconnected vessels, and the owthrough any segment depends not only on the ow resistance of that segmentbut also on the resistance of other vessels connected to it in series and inparallel.

Theoretical modeling, in combination with experimental studies, has thepotential to synthesize several observed or hypothesized mechanisms into aunied mathematical framework. The model can then be used to predict theoverall behavior of the system, taking into account the interactions betweendierent mechanisms occurring at the level of individual cells or segmentsand the interactions that arise in a network of interconnected segments. Inthis chapter we rst introduce the Hagen-Poiseuille model for ow throughducts, which provides a reasonable estimation for blood ow through vessels.Then we present two models for the distribution of hemodynamic parametersin the human retina.

2.1 Flow in the arteriolar tree

Flow in the arteriolar tree is supposed to be of Hagen-Poiseuille type.Poiseuille ows are generated by pressure gradients, with application pri-marily to ducts. They are named after J.L.M. Poiseuille (1840), a Frenchphysician who experimented with low-speed ow in tubes.Consider a straight duct of arbitrary but constant shape. There will be anentrance eect, i.e. a thin initial shear layer and core acceleration (see Fig-ure 2.1). The shear layers grow and meet, and the core disappears within afairly short entrance length,Le.For x > Le the velocity becomes purely axial and varies only with the lat-eral coordinates, so v = w = 0 and u = u(y, z) (see Figure 2.2). The ow isthen said to be fully developed. For fully developed ow the continuity and

23

Figure 2.1: Flow in the entry region of a tube. (Figure taken from [22])

momentum equations for incompressible ows reduce to:

∂u

∂x= 0

0 = −∂p∂x

+ µ

(∂2u

∂y2+∂2u

∂z2

)0 = −∂p

∂y= −∂p

∂z

These equations indicate that the total pressure p is a function only of x forthis fully developed ow. Further, since u does not vary with x, it followsfrom the x-momentum equation that the gradient ∂p/∂x must only be anegative constant. The basic equation of fully developed duct ow is thus:(

∂2u

∂y2+∂2u

∂z2

)=

1

µ

∂p

∂x= const (2.2)

Note that the acceleration terms vanish here, since the ow is very slow.

24

Figure 2.2: Fully developed duct ow. (Figure taken from [22])

Flow through a circular pipe

The ow through a circular pipe with radius R was rst studied by Hagen(1839) and Poiseuille (1840). The Laplace operator in polar coordinatesunder the hypothesis of radial symmetry and axial invariance reduces to:

∇2 =1

r

d

dr

(rd

dr

)and, under these hypothesis, the solution of the fully developed equationow, Eq. (2.2) is

u(r) =1

4µ(∂p

∂x)r2 + C1lnr + C2.

Since the velocity cannot be innite at the centerline we reject the logarithmterm and set C1 = 0. The no-slip condition is satised by setting C2 = 1

4 .The pipe-ow solution is thus:

u = −dp/dx4µ

(R2 − r2),

so that the velocity distribution in fully developed laminar pipe ow is aparaboloid of revolution about the centerline (Poiseuille paraboloid, see Fig-ure 2.3).

The total volume rate of ow Q is

Q =

∫section

udA,

which for the circular pipe gives

Qpipe =πR4

8µ

(−dpdx

). (2.3)

25

Figure 2.3: Parabolic ow in a circular pipe.

The mean velocity is dened by v = Q/A and gives, in this case

v =QpipeπR2

.

Finally, the wall shear stress is constant and given by

τw = µ

(−dudr

)w

=1

2R

(−dpdx

)=

4µv

R. (2.4)

2.1.1 Model equations

It is assumed that blood ow conforms to Hagen-Poiseuille's law in each ves-sel channel through consecutive bifurcations of the retinal microvasculature,and that the movement of material across the exchange vessels is balancedbetween blood and tissue.Hagen-Poiseuille's law indicates that the decrease in pressure ∆P againstow Q(r) along a branch of radius r and length L(r) can be written as:

∆P =8µ(r)L(r)Q(r)

πr4, (2.5)

where µ(r) is the apparent viscosity of blood that depends on the size ofthe vessel, and is supposed to follow a mathematical expression proposed byHaynes in [9],

µ(r) =µ∞

(1 + δ/r)2, (2.6)

where µ∞ is the asymptotic blood viscosity, set to 3.2·102 Poise and δ = 4.29.Blood ow exerts a tangential force that acts on the luminal surface of theblood vessel as τw(r) wall shear stress

τw(r) = µ(r)γw(r), (2.7)

γw(r) =4Q

πr3=

4v

r(2.8)

26

where γw(r) is the shear rate at the wall surface.For the conservation of ow, as represented in Figure 2.5, for each bifurcationnode the inow must be the same as the outow:

Qin = Qout (2.9)

where Qin is the total inow and Qout is the total outow.We observe that the vascular tree has analogies with a classical electricalcircuit: we can interpret the blood ow Q through vessels as the intensity ofcurrent I, the pressure drop ∆P as the potential dierence ∆V , and nallythe conductance of the vessel as the conductance of the circuit, that is theinverse of the resistance R.

Figure 2.4: Electrical circuit which represents an dichotomous arterial vas-culature.

Figure 2.5: Conservation of ow (i.e. intensity of current) in a generic net-work tree: in each node i the outow in the adjacent nodes must equal theinow so that we have

∑j∈Adj(i) Iij = Ii

27

Remembering the Hagen-Poiseuille pressure drop equation:

∆P =8µL

πr4Q ,

and comparing it with the equation governing the current ow through elec-trical network

∆V = Ri

we can see that the expression 8µLπr4 is an equivalent of the resistance R for

blood ow. We dene its inverse:

G =πr4

8µL(2.10)

as the conductance of the network. To be precise, in each bifurcation node i,where a vessel ends splitting into two daughter branches ij directed to nodesj, the vessel conductance (vessel connecting node i and node j) is dened asfollows:

Gij =πr4

ij

8µijLij(2.11)

where rij is the radius of the vessel at generation i, Lij the vessel length (thelengths of the two branches are equal for symmetry), and µij the viscosityof the vessel.Extending these node properties to the whole tree level and combining itwith the conservation of ow we get an equivalent of Kirchho law for bloodow: ∑

j∈Adj(i)

Qij =∑

j∈Adj(i)

Gij(Pi − Pj) = 0, (2.12)

where Adj(i) is the set of indexes of network nodes adjacent to node i, Piand Pj are the pressures at the nodes i and j respectively. Subsequently, wecan compute the blood ow in the single vessels Qij using:

Pi − Pj = ∆Pij =

(πr4

ij

8muijLij

)−1

Qij

=⇒ Qij = ∆PijGij . (2.13)

Boundary nodal pressures are required to start the computation. Theboundary inlet node is the artery of generation 1, where the blood owenters the network. The boundary outlet nodes are the nodes where theblood ow exits the network, in our case the capillaries.

In the case of binary network at each bifurcation nodes, because of thesymmetry of the network, the blood ow divides itself equally into the twodaughter branches:

28

Q1(r1) = 2Q2(r2) = ... = 2g−1Qg(rg),

and

vg = 2−(g−1)

(r1

rg

)2

v1,

where r1 and v1 are the radius and mean ow velocity of the trunk vessel ofgeneration 1 and g is the generic g − th ospring.

29

2.1.2 Numerical results

Using the model equations illustrated in the previous sections, we have com-puted the values of hemodynamic parameters using the binary tree and DLAnetwork, respectively. The blood pressure at the rst artery was estimatedby considering the hydrostatic and frictional pressure losses from the aortato the central retina artery, and xed at a value of 40mmHg. We considertwo conditions for the blood pressure at the outlet of the system:

• in a rst condition, we assume that 60% of nal branches goes into deepcapillary layer and that 30% goes into supercial capillary layer. Weincrease the length of these branches by a length equal to the distancebetween the plane of the arterial system and the respective planes ofthe capillaries. We set the partial pressure equal to 20mmHg in allnal branches.

• in a second condition, we set, in each nal branch, a xed blood pres-sure at the outlet of the system equal to 20mmHg.

In the binary model, we use both cases for the pressure at the outlet of thesystem. In the DLA model, we use only the rst condition in its center andwe set the branches in its periphery as branches ending at the level of thearterial system. If we compare the arterial tree to an electric circuit, wetreat the coupling between the capillary bed and arterial tree by insertingin the branches of the system that descend to this an equivalent resistancein parallel, which is calculated by comparing the capillary bed to a circuitformed by a 600 square mesh. The calculation of the equivalent resistance isdescribed in the appendix. In this study we consider an equivalent resistanceequal to 1e11g/cm4s, which we have obtained by considering the capillariesof radius 2.5µm and length 50µm. Moreover, the coupling between arterialtree and surrounding tissue in the branches that end in the arterial layeris treated inserting an equivalent resistance. We set this value equal to5e8g/cm4s.

The results are shown in the next gures. We compare the values ofmean velocity and mean ow rate that we get in the vascular system withdata measured in Riva et al [15].

30

Figure 2.6: Distribution of mean blood pressure in the binary tree modelas a function of vessel diameter. These values are obtained by imposing anoutlet pressure equal to 20 mmHg .

Figure 2.7: Distribution of mean blood pressure in the binary tree model as afunction of vessel diameter. These values are obtained by imposing an outletpressure equal to 20 mmHg and requiring that 60 % of nal branches goesinto deep capillary layer, 30 % goes into the supercial capillary bed the otherpart terminates in the arterial layer. The uctuations of mean pressure incorrespondance with each value of the vessel diameter are caused by increasedresistance of some nal branches to make the coupling to capillary beds.

31

Figure 2.8: Distribution of mean blood pressure in the DLA system as afunction of vessel diameter. These values are obtained by imposing an outletpressure equal to 20 mmHg and requiring that peripherical nal branchesterminate in the arterial layer, while internal nal branches are such that60 % of these goes into deep capillary bed, 30 % goes into the supercialcapillary bed. The high uctuations of mean pressure around the diametervalue of 30µm are caused by the existence of nal branches near the inputof DLA system.

Figure 2.9: Distribution of mean velocity in the binary tree model comparedwith experimental data in Riva et al. [15] as a function of vessel diameter .These values are obtained by imposing an outlet pressure equal to 20 mmHg.

32

Figure 2.10: Distribution of mean velocity in the binary tree model comparedwith experimental data in Riva et al. [15] as a function of vessel diameter .These values are obtained by imposing an outlet pressure equal to 20 mmHgand requiring that 60 % of nal branches goes into deep capillary layer, 30% goes into the supercial capillary bed and the other part terminates inthe arterial layer. The uctuations of mean velocity in correspondance witheach value of the vessel diameter are caused by increased resistance of somenal branches to make the coupling to capillaries beds. The highest valuesof mean velocity correspond to branches that end in the arterial layer.

33

Figure 2.11: Distribution of mean velocity in the DLA model comparedwith experimental data in Riva et al. [15] as a function of vessel diameter.These values are obtained by imposing an outlet pressure equal to 20 mmHgand requiring that peripherical nal branches terminate in the arterial layer,while internal nal branches are such that 60 % of these goes into deepcapillary bed, 30 % goes into the supercial capillary bed. The highestvalues of mean velocity are located in the central branches that end in thearterial layer. The results that are obtained with the DLA model are incloser agreement with the experimental data than those obtained with thebinary tree model.

34

Figure 2.12: Distribution of mean ow rate in the binary tree model com-pared with experimental data in Riva et al. [15] as a function of vesseldiameter. These values are obtained by imposing an outlet pressure equalto 20 mmHg .

Figure 2.13: Distribution of mean ow rate in the binary tree model com-pared with experimental data in Riva et al. [15] as a function of vesseldiameter. These values are obtained by imposing an outlet pressure equal to20 mmHg and requiring that 60 % of nal branches goes into deep capillarylayer, 30 % goes into the supercial capillary bed the other part terminatesin the arterial layer. The uctuations of mean ow rate in correspondancewith each value of the vessel diameter are caused by increased resistance ofsome nal branches to make the coupling to capillaries beds.

35

Figure 2.14: Distribution of mean ow rate in the DLA model comparedwith experimental data in Riva et al. [15] as a function of vessel diameter.These values are obtained by imposing an outlet pressure equal to 20 mmHgand requiring that peripherical nal branches terminate in the arterial layer,while internal nal branches are such that 60 % of these goes into deepcapillary bed, 30 % goes into the supercial capillary bed. The highestvalues of mean ow rate are located in the central branches that end in thearterial layer. The results that are obtained with the DLA model are incloser accordance with the experimental data than those obtained with thebinary tree model.

36

The wall shear stress at the precapillary vessels increases, since the ap-parent viscosity increases due to the geometrical obstacle encountered bythe red blood cells owing in these narrow channels. The wall shear stressof the vessels at the pre-equator and equator region is signicantly higherthan that at the periphery region, see [6]. This is reasonable because theuid is owing outward from the center of the retina, hence the further fromthe center the lower the pressure is. As a result, the driving force for theow in the arteriolar branches at the pre-equator region is higher than atthe periphery, hence the wall shear stress is higher.This leads to an observation on the relationship between the high wall shearstress and the vessel wall thickness of arterial vessels near the pre-equatorand equator regions.The wall of retinal arteries near the optic disc (pre-equator region) comprisesve to seven layers of smooth muscles. At the equator and periphery, how-ever, the arterial wall has only two or three and one or two muscle layers,respectively. This seems to suggest that the vessels at the pre-equator andequator regions have adapted themselves by increasing their wall thickness(i.e., smooth muscles) to sustain the higher wall shear stress.

Figure 2.15: Distribution of mean shear stress in the binary tree system asa function of vessel diameter. These values are obtained by imposing outletpressure equal to 20 mmHg .

37

Figure 2.16: Distribution of mean shear stress in the binary tree system asa function of vessel diameter. These values are obtained by imposing outletpressure equal to 20 mmHg and requiring that 60 % of nal branches goesinto deep capillary layer, 30 % goes into the supercial capillary bed theother part terminates in the arterial layer. The uctuations of mean shearstress in correspondance with each value of the vessel diameter are caused byincreased resistance of some nal branches to make the coupling to capillarybeds.

38

Figure 2.17: Distribution of mean shear stress in DLA system as a functionof vessel diameter. These values are obtained by imposing an outlet pressureequal to 20 mmHg and requiring that peripherical nal branches terminatein the arterial layer, while internal nal branches are such that 60 % of thesegoes into deep capillary bed, 30 % goes into the supercial capillary bed.The highest values of mean shear stress are located in the central branchesthat end in the arterial layer.

39

Figure 2.18: The plot represents the spatial distribution of mean blood pres-sure in the DLA system. The colors represent the values of this. The pres-sure decreases along the main artery from the center to the periphery, toreach 20mmHg. The lowest pressures at the center of system are caused byimposing the Dirichlet boundary conditions.

40

Figure 2.19: The plot represents the spatial distribution of mean velocity inthe DLA system. The colors represent the values of this. The mean velocitydecreases along the main artery reaching very low levels in the periphery.The highest values of the velocity correspond to the vessels that end inthe arterial plane and this fact is due to the high resistance applied to thebranches that descend into the capillary beds.

41

Figure 2.20: The plot represents the spatial distribution of mean ow ratein the DLA system. The colors represent the values of this. The mean owrate decreases along the main artery reaching low values in the periphery.

42

Figure 2.21: The plot represents the spatial distribution of mean shear stressin the DLA system. The colors represent the values of this. The highest valesof mean shear stress correspond to the vessels in which there are the highestvalues of mean velocity. Also, between these vessels, the value of mean shearstress is greatest if the radius is smallest.

43

2.2 Flow in Capillary Plexi and in Interstitial Tis-

sue

The retinal capillaries are embedded in the retinal tissue, but these domainsare very dierent from each other. To describe these, we use the double-continuum approach that we presented in Ch.1 in which they are treatedas two separate continua coupled by exchange functions. In this section webriey describe the equations that govern the uid ow in these domains.

Retinal Tissue

Since we consider interstitial tissue as a porous medium, the ow velocity ofthe interstitial uid can be described by Darcy's law:

−→v F = −Kµ

(∇P ), (2.14)

where −→v F is the Darcy velocity, K is the intrinsic permeability tensor andµ the uid viscosity of the uid phase.The conservation of mass is expressed by the following equation:

∇ · (ρT−→v F ) + q − qF = 0 (2.15)

where ρT is the apparent mass density, i.e. the intrinsic mass density of uidper unit volume (ρtissue ) multiplied by the volume fraction of tissue withinthe model domain (fT = Vtissue/Vtot), q represents the external sources orsink (e.g. the inuence of the lymphatic system, that can be seen as a sink),qF is the coupling variable for the ow between the two continua. Sincewe consider a incompressible uid phase and a constant tissue porosity φ,the temporal variation of the product of φ and ρT is not considered in thecontinuity equation.

Retinal Capillary

As in the case of the retinal tissue, the capillary bed can be treated asa porous medium, so that Darcy's law may be applied to determine theblood ow velocity. The ow of blood can be described with the followingcontinuity equation:

∇ · (ρC−→v F ) + qF = 0 (2.16)

where ρC is the apparent mass density, i.e. the intrinsic mass density of uidper unit volume (ρblood) multiplied by the volume fraction of tissue within themodel domain (fC = Vcapillary/Vtot), qF is the coupling variable for the owbetween the two continua. Since we consider a incompressible uid phaseand a constant tissue porosity φ = 1, the temporal variation of the productof φ and ρC is not considered in the continuity equation.

44

Coupling Tissue-Capillary

The exchange of uid across the capillary walls into the retinal tissue andvice-versa, is governed by an exchange term, the transfer or coupling func-tion qF , that depends on the uid pressure gradient across the vessel wall.According to Starling's law, net uid ow across a vessel wall is given by:

qF = ρmolLpAvesselVtissue

(Pc − Pis) (2.17)

where Lp is the hydraulic conductivity of the vessel wall, Avessel/Vtissue is thesurface area of the retinal capillaries per unit volume of tissue and ρmol is theintrinsic mass density of uid that passes from capillary into the interstitiumacross the microvascular wall. The term within parenthesis is called thetransmural pressure, where pc and pis are the uid pressure in capillariesand interstitial space, respectively.

2.2.1 Calculation of Permeability in Capillary Bed

In this section, we present the computation of the continuum intrinsic perme-ability tensor in the capillary bed, that depends mainly on the connectivityof the capillary segments and their size. We divide the capillary domaininto 600 rectangular cuboid subvolumes (see Fig.2.22 ) and we calculate thepermeability tensor as in Reichold et al.[14]. In the x-direction (in a similarway in the y-direction) we have:

Kx =Fx,REV µLx,REV

| PC − PA | Ay,z,REV(2.18)

where Fx,REV is the mass ow between the two faces normal to the x-axis,µ is the dynamic viscosity of blood, PA and PC are the pressures in A and Crespectively, Ay,z,REV is the cross section area of REV parallel to y-z-plane.

The volume of blood is:

Vblood ' 2πr2Lx − (2r)3 (2.19)

and the porosity of blood is:

φblood =VbloodVREV

=2πr2Lx − (2r)3

L2xδ

. (2.20)

To calculate the value of Kx, we set an arbitrary pressure boundarycondition PA and PC and no-ow boundary conditions on the face parallelto x-z-plane and x-y-plane . Since we suppose a linear trend for the pressurebetween A and C, the pressure value in B is (PC + PA)/2. because ofsymmetry, we consider the cylinder between B and C and assuming Hagen-Poiseuille blood ow, Fx,REV is:

Fx,REV =πr44P4µLx

(2.21)

45

Figure 2.22: The REV. We consider Lx=50µm, δ=30µm and r=2.5µm.

where

4P =| PB − PC |=| PC − PA |

2.

Hence the permeability in the subvolume in x-direction is:

Kx =πr4

8δLx. (2.22)

The numerical value of Kx is 0.01 µm2.

46

Chapter 3

Models for Solute Transport

and Delivery

In this Chapter we present the models which describe the transport of asolute from blood ow through the retina tissue. First, we describe solutetransport across the arterial wall using a Wall-Free model and a Multilayermodel, then we present a model for solute transport in the capillary bedand in the retinal tissue. At the end of each section we adapt the modelequations in the case the solute is constituted by oxygen.

3.1 Solute transport across the vessel wall

A schematic radial cross-section of a vessel consists of: 1) a red blood cell-richcore (RBC), 2) a red blood cell-depleted plasma layer (PL), 3) endothelialvascular wall (ET), 4) smooth muscle layer (SMC) and 5) tissue space. Thelumen comprises the RBC and PL layers (see Fig. 3.1).In our study we consider that the lumen of a vessel is not divided actually inRBC and PL, but is a unique continuum. Each blood vessel is considered as acylinder of given radius Rv and length Lv in a three-dimensional system withcylindrical coordinates. The z axis corresponds to the axis of the cylinder,as in Fig. 3.2. The general geometrical structure of the model consists inconcentric cylinders as shown in Fig .3.2.

We denote by ΩL the lumen, ΩET the endothelial wall, ΩSMC the smoothmuscle layer and ΩT the tissue layer. We indicate by Γi the part of ∂Ωi

corresponding to the common interface between two distinct domains, by−→n i the unit outward normal vector with respect to ∂Ωi and Ri the distanceof the interface Γi from the z axis. Moreover, in the following we will denoteby ΓL,in and ΓL,out the inlet and outlet section of the lumen, respectively.

47

Figure 3.1: Schematic cross-section of the model geometry showing the multi-layer structure of the vessel: RBC-rich core (CR), RBC-depleted plasmalayer (PL), endothelial vessel wall (ET), smooth muscle cell (SMC), andtissue space (TS).

Figure 3.2: The considered subdomains and the partitioning of the bound-aries for the arterial wall model.

48

3.1.1 The Wall-Free Model

In this section we consider only a lumen as a physical domain ΩL, while thepresence of the arterial wall is expressed by a "transfer" boundary condition.We suppose that the solute, in addition to being transported and diusedalong the vessel, also lters throughout the wall of the vessel itself.

The equations that describe the transport and diusion of solute in ΩL

are:

−D[

1

r

∂

∂r

(r∂P

∂r

)+∂2P

∂z2

]+

∂

∂z(vzP ) = 0 in ΩL, (3.1a)

−D∂P∂r· −→n L = γ(P − Pwall) on ΓL , (3.1b)

P = Pinlet on ΓL,in , (3.1c)

J · −→n L,out − v+n P = αP − Jout · −→n L,out on ΓL,out. (3.1d)

whereD is the diusion coecient of the solute that we suppose constant,A the constant cross-section area of ΩL,

−→v the velocity of the blood, P isthe solute partial pressure, Pwall is the solute partial pressure in the vascularwall, γ is the solute permeability across the lumen-wall interface, J is a uxof density mass and ∂Ω represent a wall of the vessel. We consider a Dirichletcondition at the start of vessel (ΓL,in ) and a Robin condition at the end ofthe vessel (ΓL,out ) [10]: Pinlet is the partial pressure at the inlet of domain,Jout is the ux of density mass at the distal of domain and α is a constantvelocity .We dene

vn = −→v · −→nv+n = (vn + |vn|)/2v−n = (vn − |vn|)/2vn = v+

n + v−n .

It is convenient to dene a new variable Pav = Pav(z), the average of Pon the section area:

Pav =

∫ 2π0

∫ R0 rP (r, z)drdθ

πR2. (3.2)

We integrate Eq.(3.1a) on the section, obtaining:

49

∫ 2π

0

(∫ R

0−D ∂

∂r

(r∂P

∂r

)dr +

∫ R

0−D∂

2P

∂z2rdr +

∫ R

0

∂

∂z(vzP ) rdr

)dθ = 0.

Using the fundamental theorem of calculus in the rst term, we obtain:

−DR ∂P

∂r

∣∣∣∣ΓL

+

∫ R

0−D∂

2P

∂z2rdr +

∫ R

0

∂

∂z(vzP ) rdr = 0,

and using the boundary conditions, Eq.(3.1b),

Rγ(P − Pwall) +

∫ R

0−D∂

2P

∂z2rdr +

∫ R

0

∂

∂z(vzP ) rdr = 0. (3.3)

Using Eq. 3.2 and assuming that the variation of v is small along z in theblood vessel, we replace vz with vmean, and we approximate P on ΓL withPav. We obtain the following averaged model for the transport of solute inthe vessel:

∂J

∂z+

2

R2Rγ(Pav − Pwall) = 0 in ΩL, (3.4a)

J = −D∂Pav∂z

+ vmeanPav in ΩL, (3.4b)

Pav = Pinlet on ΓL,in, (3.4c)

J · −→n L,out − v+n Pav = αPav − Jout · −→n L,out on ΓL,out, (3.4d)

To describe the solute transport in each vessel of the arterial tree we useEq. (3.4a). Eqs (3.4c)-(3.4d) are used as boundary conditions at the inletand outlet of the arterial tree, respectively. At each bifurcation node s, weassume the continuity of mass uxes through each section:

AiJi · −→n i = AjJj · −→n j +AkJk · −→n k, (3.5)

where i,j and k are the names of vessels that have s as common node, Arepresents the section area, J represents the density mass ux and −→n theunit outward normal vector with respect to boundary of the vessel at thebifurcation node.We solve numerically the convection-diusion-reaction equations of the Free-Wall Model using the FEM and the Scharfetter-Gummel method. Thesemethods are described in the appendix.

50

3.1.2 The Multilayer Model

In this section we study the same problem as in the previous section, butwhit a geometry that is more similar to the real one. We consider, in thismodel, a lumen, an endothelium layer and a smooth muscle cell layer. Thepresence of the tissue layer is expressed by a boundary condition (see Fig.3.2)We suppose that the solute transport occurs within the arteriole space viaconvection and diusion through the plasma. The general mass balance forthe solute inside the arteriole is given by:

−2π(rDL∂PL∂r

)|ΓL+∂I(Pav,L)

∂z= 0, in ΩL

where

I(Pav,L) = πR2L(−DL

∂Pav,L∂z

+ vmeanPav,L),

Pav,L =

∫ 2π0

∫ RL

0 PL(r, z)rdrdθ

πR2L

,

and where vmean is the mean of blood velocity in the lumen, as in the pre-vious section.In this model, the transport of solute in the arterial wall is controlled by dif-fusion through this region and consumption by endothelial cells and smoothmuscle cells. The rate of solute consumption follows the Michaelis-Mentenkinetics. So the transport in the wall is described by the following equation:

−1

r

∂

∂r

(rDi

∂Pi∂r

)+

Km,iPiK1/2,i + Pi

= 0 in Ωi (3.6)

where Pi = Pi(r) is the radial partial pressure of solute in the arteriole wall,Di the diusion coecient, Km,i and K1/2,i the constants of the Michaelis-Menten law. The subscript i indicates whether we are in endothelium layer(ET) or smooth muscle cells layer (SMC). We assume continuity of the solutepartial pressure and the mass ux at each boundary interface, so that theequations that describe the transport and diusion of solute in lumen andwall are:

51

−2π(rDL∂PL∂r

)|ΓL+∂I(Pav,L)

∂z= 0 in ΩL, (3.7a)

−1

r

∂

∂r

(rDE

∂PE∂r

)+

Km,ETPETK1/2,ET + PET

= 0 in ΩET , (3.7b)

−1

r

∂

∂r

(rDSMC

∂PSMC

∂r

)+

Km,SMCPSMC

K1/2,SMC + PSMC= 0 in ΩSMC , (3.7c)

Pav,L = Pinlet on ΓL,in, (3.7d)

J · −→n L,out − v+n Pav,L = αPav,L − Jout · −→n L,out on ΓL,out, (3.7e)

−DL∂PL∂r

+DET∂PET∂r

= 0 on ΓL, (3.7f)

PET = Pav,L(z = Lv/2) on ΓL, (3.7g)

−DET∂PET∂r

+DSMC∂PSMC

∂r= 0 on ΓET , (3.7h)

PET = PSMC on ΓE , (3.7i)

DSMC∂PSMC

∂r+ hWPSMC = hWPT on ΓSMC . (3.7j)

We use the FEM and Scharfetter-Gummel method with the assumptionof continuity of mass ux through each vessel section to solve the systemequations in (3.7) corresponding to the lumen. The FEM and Scharfetter-Gummel method are described in the Appendix.The equations that describe the solute transport in the endothelium andsmooth muscle cells layers, Eqs. (3.7b)-(3.7c), are non-linear, for this reasonwe use a xed point method to linearize and solve them. We approximatethe variable Pi with its mean across the radius of the layer, i.e.:

PET ' Pav,ET :=1

RET −RL

∫ RET

RL

PET (r)dr in ΩET ,

PSMC ' Pav,SMC :=1

RSMC −RET

∫ RSMC

RET

PSMC(r)dr in ΩSMC ,

52

and we obtain:

−1

r

∂

∂r

(rDET

∂PET∂r

)+ γ(Pav,ET )PET = 0, in ΩET , (3.8)

−1

r

∂

∂r

(rDSMC

∂PSMC

∂r

)+ γ(Pav,SMC)PSMC = 0, in ΩSMC . (3.9)

where

γ(Pav,i) =Km,i

K1/2,i + Pav,iwith i = ET, SMC

Eqs (3.8) and (3.9) have the form of the so-called modied Bessel equationsof order 0 (Kelvin's equations of order 0). Their solutions are [1]:

PET = AI0

(r/

√DET

γ(Pav,ET )

)+BK0

(r/

√DET

γ(Pav,ET )

)in ΩET ,

PSMC = CI0

(r/

√DSMC

γ(Pav,SMC)

)+DK0

(r/

√DSMC

γ(Pav,SMC)

)in ΩSMC ,

where I0 and K0 are the modied Bessel functions of the rst and secondkind of order 0, and A, B, C and D are integration constants.

The Multilayer model can be seen as a coupling model which combines asolute transport model in the vessel lumen (Eq.(3.7a)) and a solute transportmodel in the vessel wall, ET and SMC (Eqs.(3.7b) and (3.7c)). The inputdata for the model are: the inlet and outlet boundary conditions for thelumen, the mean partial pressure in ET and SMC and the value of thederivative of partial pressure in lumen in r-direction. We solve the lumenmodel nding all its parameters. These are used to communicate with thewall model until, after performing a xed point iteration over the iterationcounter k, the following conditions are satised,

‖P (k+1)av,ET − P

(k)av,ET ‖∞

‖P (k)av,ET ‖∞

≤ Tol,

‖P (k+1)av,SMC − P

(k)av,SMC‖∞

‖P (k)av,SMC‖∞

≤ Tol.

Then the found values are used to communicate with the lumen model.This process is repeated until

53

‖P (t+1)av,L − P

(t)av,L‖∞

‖P (t)av,L‖∞

≤ Tol.

t represents the outer iteration counter, i.e. the number of iterations of theoverall process.

3.1.3 The case of solute oxygen: free oxygen and oxygen

carried by oxyhemoglobin

In blood there is a reversible dynamic reaction between the oxygen carrier,Hemoglobin, and free oxygen in the plasma. Hemoglobin (Hb) is a proteincontained in the red blood cells that transports oxygen from the lungs tothe peripheral tissues of the body. In this section we present a model in-cluding both the free oxygen dissolved in plasma and the oxygen carried byhemoglobin, as in [11]. The governing equation in plasma is:

∇ · (−Dp∇[(1−HD)Cp] +−→V p(1−HD)Cp) = J (3.10)

where Cp is the oxygen concentration in plasma, HD is the hematocrit (thevolume percentage (%) of RBCs in blood), Dp is the diusivity of free oxygen

in plasma,−→V p is the velocity of plasma and J is the released oxygen ux by

the RBCs. (1−HD)Cp represents the percentage of the free oxygen.It is assumed that free oxygen in the RBCs can easily pass through themembrane of erythrocyte and exchange with outside plasma, while oxygenbound to hemoglobin only exists in RBCs.In RBCs, we have:

∇ · (−DHb∇(HDCHbS)−Dc∇(HDCc) +−→V rbcHD(Cc + CHbS)) =

= −J (3.11)

where Cc is the free oxygen in RBCs, DHb is the diusivity of oxyhemoglobin

in RBCs, Dc is the diusivity of free oxygen in RBCs,−→V rbc is the velocity

of RBCs, CHb is the oxygen-carrying capability of hemoglobin in blood andS is the oxyhemoglobin saturation function expressed by the Hill equation:

S = Pnc /(Pnc + Pn50) (3.12)

where P50 is the half partial pressure of oxygen saturation hemoglobin andn is the Hill exponent. HD(Cc+CHbS) represents the percentage of the freeoxygen and the oxygen combined with Hb.

54

According to Henry's law, the free oxygen concentration (Ci) and corre-sponding partial pressure (Pi) are related by

Ci = αiPi

where αi is the oxygen solubility coecient in the relevant uid. ThereforeEqs.(3.10) and (3.11) become:in plasma

∇ · (−Dp∇[(1−HD)αpPp] +−→V p(1−HC)αpPp) = J, (3.13)

in RBCs:

∇ · (−DHb∇(HDCHbS)−Dc∇(HDαcPc) +−→V rbcHD(αcPc + CHbS)) =

= −J. (3.14)

Adding Eq.(3.13) and Eq.(3.14), we obtain:

∇ · (−→V rbcHD(αcPc + CHbS) +

−→V p(1−HC)Cp) =

= ∇ · (DHb∇(HDCHbS) +Dc∇(HDαcPc) +Dp∇[(1−HD)αpPp]). (3.15)

Assuming that there is no relative movement between plasma and RBCs,i.e., −→

V p =−→V rbc =

−→V b

and that the diusivity of the free oxygen in RBCs is the same as that inplasma,Dc = Dp and that Pp = Pc = Pb and αp(1−HD) + αcHD = αb, wenally get:

∇ · (−→V p(Pb +

HDCHbαb

S) = ∇ · (Dp∇Pb +DHbHDCHb

αb

∂S

∂Pb∇Pb)

and

∇ · [−(Dp + DHbHDCHbαb

∂S∂Pb

)∇Pb) +−→V p(Pb + HDCHb

αbS)] = 0,

∇ · [−(Dp + DHbHDCHbαb

∂S∂Pb

)∇Pb +−→V pPb(1 + HDCHb

αb

Pn−1b

(Pnb +Pn

50))] = 0.

We let:

η(Pb) = Dp +DHbHDCHb

αb

∂S

∂Pb, ξ(Pb) = 1 +

HDCHbαb

Pn−1b

(Pnb + Pn50)

55

and obtain:∇ · [−η(Pb)∇Pb +

−→V p(Pbξ(Pb)] = 0. (3.16)

Since (3.16) is a non-linear dierential equation, we use a xed pointmethod to linearize and solve it. We iterate the process until the relativeerror is smaller than a given tolerance Tol, i.e.:

‖P (k+1)b − P (k)

b ‖‖P (k)

b ‖≤ Tol.

Using the Wall-Free Model, we consider this boundary condition on theinterface ΓL:

−−→n · (η(Pb)∇Pb) = P0(Pb − Pwαwαb

) (3.17)

where −→n w is the normal vector of Γw, P0 denotes the oxygen permeabil-ity across the lumen-wall interface, Pw the oxygen partial pressure in thevascular wall and αw the oxygen solubility coecient in the vessel wall.

The values of the parameters used in the previous equations are given inTab. 4.1

56

3.2 Solute Transport in Capillary Beds and in In-

terstitial Tissue

In this section, we want to propose a method to treat the transport of thesolute in the capillary beds and in the interstitial tissue. As discussed inCh.1, we use a double continuum approach to describe solute transport inthe capillary layer.The transport for the solute in the retinal tissue is described by the followingequation:

∂(CT fTφ)

∂t+∇ · (CT fT−→v F − fTDeff∇CT ) + CT − qT = 0 (3.18)

where CT is the solute concentration in the tissue, fT is the tissue volumefraction, φ is the porosity, −→v f is the Darcy velocity, Deff is the eectivediusion coecient of the solute, rT represents the external sources or sinkand qT is the transport coupling variable. The second term describes theadvection and diusive transport of the solute within the tissue.The transport for the solute in the retinal capillary bed is described by thefollowing equation:

∂(CCfCφ)

∂t+∇ · (CCfC−→v F − fCDeff∇CC) + rC + qT = 0 (3.19)

where CC is the solute concentration in capillary, fC the capillary volumefraction, φ is the porosity, −→v f is the Darcy velocity, Deff is the eectivediusion coecient of the solute, rC represents the external sources or sinkand qT is the transport coupling variable. The second term describes the ad-vection and diusive transport of the solute within the tissue. In both cases,the eective diusion coecient is dened as the product of the tortuosityfactor, the porosity of the medium and the diusion coecient of solute inthe uid phase.

Solute transport between the tissue and capillary continuum is describedby the coupling function qF , based on the characteristics of the transportacross the capillary wall that depends on the relative concentration gradient.The Starling equation describes the advective and diusive transport of thesolute across the vessel wall:

qT = PAvesselVtissue

(Cc − CT ) (3.20)

where P is permeability of the capillary wall, Avessel/Vtissue is the surfacearea of the retinal capillaries per unit volume of tissue.

57

3.3 Solute transport in the tissue

In this section, we study the solute partial pressure distribution in the retinaltissue. The model presented here is based on those proposed in [3] and in[4], in which it is supposed that the retina is divided into eight layers, asdescribed in Ch.1, each with a distinct solute consumption or supply rate.We suppose that, in the retinal tissue, the transport of solute is only due todiusion. Therefore, in steady state conditions, the change in solute partialpressure P is given by the equation:

∇ · J − q + s = 0,

J = −D∇P,(3.21)

where D is the solute diusion coecient, q is the solute consumptionterm and s is a delivery term.According to the Michaelis-Menten equation, the consumption term for so-lute q is given by:

q =PKm

P +K1/2, (3.22)

where Km is the maximal rate of solute consumption and K1/2 is the partialpressure of solute at half maximal consumption speed. The amount of solutetransported from blood to tissue, according to the Fick principle, is givenby:

s = Q(αPblood − β(P )P ), (3.23)

where Q is the blood ow rate, Pblood is the partial pressure of solute inarterial blood, α is a constant and β(P ) is a function of P that describes thesource term of oxygen in tissue.

We assume to neglect solute diusion except along the Choroid-Vitreousdirection (see Fig. 1.1). Therefore the solute diusion model (3.21) becomesa 1D model on the domain Ω = ∪j=1,..,8Ωj (see Fig. 3.3), where z is thedistance from the Choroid, Ωj are the layers introduced in Ch.1, ΓC theinterface between retina and Choroid and ΓV the interface between retinaltissue and vitreous.

Thus, we have eight dierent sub-models of the form Eq. (3.21), one foreach domain Ωj . We enforce that both solute partial pressure and ux J arecontinuous on Γj . The coecient diusion Dj is set equal in each layer andwe consider a Dirichlet condition on ΓC and a Neumann condition on ΓV .Then we obtain:

58

Figure 3.3: The domain for the model of solute transport in the tissue.

−D∂2Pj∂x2

+ qj − sj = 0 in Ωj , j = 1, .., 8,

−D∂Pj∂x

= −D∂Pj+1

∂xin Γj , j = 1, .., 7,

Pj = Pj+1 in Γj , j = 1, .., 7,

P1 = PC in ΓC ,

−D∂P8

∂x= 0 in ΓV ,

(3.24)

where Pj = P|Ωj, qj is the consumption term as in Eq. (3.22), sj is the

delivery term as in Eq. (3.23), PC the solute partial pressure at ΓC .

59

3.3.1 The O2 model for Retinal Tissue

Now we specialize the previous model using oxygen like solute. Layers 1, 3,5, and 7 are assumed to have negligible oxygen consumption. This is basedon the known properties of layer 1 and 3 from the outer retina studies andobservations in the rat inner retina [4]. Since oxygen supply and consumptionare intermingled in layer 4 and 8, the absolute level of oxygen consumptioncannot be quantied. Oxygen supply and oxygen consumption in layer 4and 8 can counterbalance each other. Therefore oxygen consumption andoxygen supply in layer 4 and 8 can be assumed to be negligible. The oxygensource term only applies to the inner retina since the outer retina does nothave any blood ow.We assume that the delivery term is of the form [16]

sj =bf60

[(Pblood − Pj) +

(Pnblood

Pnblood + Pn50

− (Pj)n

(Pj)n + Pn50

)Hbδ

α1

](3.25)

where the rst term is the source of oxygen coming from the free oxygen inthe blood and the second term the oxygen source coming from the Hb in theRBCs. Moreover, Pblood is the partial pressure in arterial blood, Hb is thehemoglobin concentration in blood, α1 is the solubility of oxygen in blood,δ is the oxygen carrying capacity of hemoglobin and P50 is the half partialpressure of oxygen saturation in Hb. The description of the hemoglobinsaturation curve (the expression within parenthesis multiplied by Hb) is thewell-known Hill equation, where n is the Hill coecient (Eq.(3.12)).

3.3.2 Numerical Solution

We solve the previous model with the FEM using a xed point procedure totreat the non-linearity of the delivery term. The value of parameters thatwe use are reported in the Table 4.1.Fig. 3.4 shows the retinal oxygen prole under normal condition (light anddark), with Pblood = PC = 80 mmHg and bf=0.2. The partial pressureis maximal at the Choroid and decreases sharply in correspondance of theouter region of inner plexiform layer (layer 6) and in correspondance of theinner segments of photoreceptors (layer 2) where there is the major oxygenconsumption.

60

Figure 3.4: Retinal oxygen proles in normal condition, where pblood was setto 80 mmHg. The yellow line represents the oxygen prole in dark conditionwhereas the blue line represents the oxygen prole in light condition.

Fig. 3.5-3.6 show the oxygen proles under two dierent hyperoxia con-ditions, with bf=0.4, Pblood = PC = 250 mmHg, 405 mmHg. We can seethat, in such condition, the oxygen consumption of level 6 has little inuencein the decrease of oxygen partial pressure in the inner retina.

Figure 3.5: Retinal oxygen proles in hyperoxia condition, where pblood wasset to 250 mmHg. The yellow line represents the oxygen prole in dark con-dition whereas the blue line represents the oxygen prole in light condition.

61

Figure 3.6: Retinal oxygen proles in extreme hyperoxia condition, wherepblood was set to 405 mmHg. The yellow line represents the oxygen prolein dark condition whereas the blue line represents the oxygen prole in lightcondition.

Fig. 3.7-3.8 show the oxygen prole that we obtain if we set all the sourceterms in the inner retina equal to 0 (avascular retina). We can see that innormal condition, the oxygen that comes from the blood vessels of the retinaplays a fundamental role in tissue oxygenation. The choroid is not able toadequately feed all the layers of the retina. Under hyperoxic condition, theoxygen from retinal vessel seems not be needed to feed retinal tissue sincethe Choroid supplies a sucient quantity of oxygen.

62

Figure 3.7: Oxygen proles in avascular retina under normal condition,where PC was set to 80 mmHg. The yellow line represents the oxygen prolein dark condition whereas the blue line represents the oxygen prole in lightcondition.

Figure 3.8: Oxygen proles in avascular retina under hyperoxia condition,where PC was set to 250 mmHg. The yellow line represents the oxygen prolein dark condition whereas the blue line represents the oxygen prole in lightcondition.

63

Chapter 4

Multiscale coupled Model

In this Chapter we present a coupled model, which combines the blood owand the solute transport models (Wall-Free and Multilayer model) with thetissue model reported in Ch. 3. We propose an iterative approach thatterminates when convergence between the models is achieved. At the end ofthis section, a summary of the model equations is presented in Fig. 4.9. Weuse both solute transport models in arterioles since we suppose that vesselswith radius under 15µm have a wall very thin and for these the layers ofSMCs and ET can be considered as a only domain. So, for these vessels,we use the Wall-Free approach. In the other vessels, we use the Multilayermodel and we suppose that for each of these vessels a relationship betweenthickness t of wall and radius r holds:

t√r = const = Tmax

√Rmax,

where Rmax is the maximum radius in arterial network and Tmax is thethickness of its wall. This value is reported in Tab.4.1.

The input data of the model are:

• the inlet and outlet blood pressure in the arterial network,

• the inlet and outlet solute concentration in the arterial network,

• the solute concentration in the vessel walls.

First, the network for blood pressure is solved and all its characteristic quan-tities are found (Mean Pressure, Mean Flow Rate, Mean Velocity and MeanShear Stress). Then, we solve the solute transport model in vessels lumenas described above and then we solve the solute transport model in wall ofvessels with radius greater than 15µm. After , we solve the solute transportin tissue model and we check if the "external" convergence is reached; if not,we iterate again the process solving the solute transport model in the arterialnetwork.

64

H

Figure 4.1: A schematic representation of the algorithm followed to implement the coupled model.

In the rst iteration of this model, we use the Wall-Free model in allvessels in order to have a initial distribution of solute values in vessels thatrespects the geometry of the arterial tree. A schematic representation of thisalgorithm is shown in Fig. 4.1.

4.1 Numerical results

In this section we show the numerical results if we consider oxygen as solute.Before doing this we describe the transfer functions that allow the couplingbetween the arterial network and the tissue model. The coupling betweensolute transport model in vessel lumen and wall has been already presentedin the Ch. 3. Moreover, we use Dirichlet conditions to calculate the oxygenpartial pressure in the vessels lumen at outlet of the arterial network unlikeRobin boundary conditions presented in the Ch. 3 in order to connect thesemodels with the model for solute transport in the tissue.The variable bf in (3.25) is dened as the blood ow that arrives in a tissueportion (assumed to be a cylinder), normalized by the weight of the portionitself. We dene Q as the harmonic mean of all mean ows of vessels that gointo the capillary beds and S as the harmonic mean of all vessel cross sections.The weight of tissue cylinder on which a vessel of section S is immersed isgiven by S ∗ tcyl ∗ ρtissue, where we set tcyl = 10−2cm and ρtissue = 1g/cm3.Therefore, bf is equal to :

bf =Q

S ∗ tcyl ∗ ρtissue. (4.1)

Moreover, we chose Pblood as the harmonic mean between all mean oxygenconcentrations of vessels that go into capillary bed. The tissue model com-municates to the arterial network the oxygen partial pressure in the vitreous(Γv in tissue domain), and the harmonic mean of oxygen partial pressure ineach capillary bed (layer 4 and 8). These values are used to dene the partialpressure of oxygen in the tissue and in the vessel wall, Pwall and PT , andto dene the boundary conditions of the vessels that go into the capillarylayers, respectively. The boundary conditions for the other vessels are setequal to PT .

Figs. 4.2-4.3-4.4 show the distribution of oxygen partial pressure in lu-men, endothelium layer and smooth muscle cells layer, respectively, of eachvessel in the arterial network. In each of these gures, the values of pressureare high in the center of the network and decrease along the four main ar-teries from center to periphery (see also Fig. 4.7). While the values of theoxygen partial pressure for the vessels with radius less then 15µm in theFig. 4.2 represent the real values of this, calculated with Wall-Free model, in

66

the other gure these values are "dummy" values set equal to zero since inthese vessels we do not consider the endothelium and the smooth muscle cellslayers. Comparing these gures, we can see that if we consider a particularvessel of the arterial tree, the oxygen partial pressure decreases passing fromthe lumen to the endothelium and nally to the smooth muscle cells layer(see also Fig.4.8).

Figure 4.2: Colours represent the value of oxygen partial pressure in vessellumen.

67

Figure 4.3: Colours represent the value of oxygen partial pressure in vesselendothelial layer.

68

Figure 4.4: Colours represent the value of oxygen partial pressure in vesselsmooth muscle cell layer.

69

Fig 4.5 shows the distribution of oxyhemoglobin saturation in vessel lu-men, dened as Eq. (3.12). This has the same trend of the distribution ofoxygen partial pressure and its values pass from about 95% to 55% from thecenter to the periphery (see Fig. 4.7).

Figure 4.5: Colours represent the value of oxyhemoglobin saturation in vessellumen.

Fig 4.6 shows the retinal oxygen prole in the tissue, in light condition,obtained by the coupling with the arterial network. The partial pressure ismaximal at the Choroid and decreases sharply in correspondence of layer 2and 6 (inner plexiform layer and inner segments of photoreceptors respec-tively). This trend is similar to oxygen prole for the retinal tissue in theliterature.

70

Figure 4.6: Retinal oxygen prole in light condition obtained with the cou-pled model. The value of Pblood represents the average of the mean partialpressure of oxygen in the vessels that reach the capillary beds and is used asinput in the tissue model.

71

Figure 4.7: The plots represent the value of oxygen partial pressure andoxyhemoglobin saturation, respectively, as a function of distance from thecenter of arterial network.

72

Figure 4.8: The plots represent the value of oxygen partial pressure in lumen and wallof three dierent vessels of the arterial tree, chosen between the vessels with radiusgreater then 15µm, as a function of distance from the axis of the vessel themself. Thedotted lines represent the interfaces between lumen-endothelium, endothelium-SMCsand SMCs-tissue, respectively. Since in arterial models we approximate the valueoxygen partial pressure in each point of the lumen by its average on the section area,we have lost the information about the oxygen prole into radial direction of the lumenand for this reason we represent it by a constant.

73

Multilayer

ModelEquations

TissueModelEquations

−2π

(rDL∂PL

∂r

)| ΓL

+∂I(P

av,L

)

∂z

=0

inΩL,

−D∂

2Pj

∂x

2+q j−s j

=0

inΩj,j

=1,..,8,

−1 r

∂ ∂r

( rDE∂PE

∂r

) +Km,ETPET

K1/2,ET

+Pav,ET

=0

inΩET,

−D∂Pj

∂x

=−D∂Pj+

1

∂x

inΓj,j

=1,..,7

−1 r

∂ ∂r

( rDSMC∂PSMC

∂r

) +Km,SMCPSMC

K1/2,SMC

+Pav,SMC

=0

inΩSMC,

Pj

=Pj+

1in

Γj,j

=1,..,7,

Pav,L

=Pinlet

onΓL,in,

P1

=PC

inΓC,

J·−→nL,out−v

+ nPav,L

=αPav,L−Jout·−→nL,out

onΓL,out,−D∂P

8

∂x

=0

inΓV.

−DL∂PL

∂r

+DET∂PET

∂r

=0

onΓL,

Wall-FreeModelEquations

PET

=Pav,L

(z=Lv/2)

onΓL,

−DET∂PET

∂r

+DSMC∂PSMC

∂r

=0

onΓET,

∂J

∂z

+2 R2Rγ

(Pav−Pwall)

=0

inΩL,

PET

=PSMC

onΓE,

J=−D∂Pav

∂z

+v m

eanPav

inΩL,

DSMC∂PSMC

∂r

+hWPSMC

=hWPT

onΓSMC.

Pav

=Pinlet

onΓL,in,

J·−→nL,out−v

+ nPav

=αPav−Jout·−→nL,out

onΓL,out.

Figure

4.9:

Summaryof

modelequationsforsolute

transport.

74

Symbol Value Description Reference

Arteriole lumenPinlet 100 mmHg Arteriole inlet O2 partial pressure /P50 26.6 mmHg Half partial pressure of O2 saturation in Hb [11]n 2.7 Hill exponent [11]CHb 0.2 mlO2/ml O2 carrying capability of Hb in blood [11]αp 2.82e-5 mlO2/ml/mmHg O2 solubility coecient in plasma [11]αc 3.38e-5 mlO2/ml/mmHg O2 solubility coecient in RBCs [11]DHb 1.5e-7 cm2/s Diusivity of Hb in RBCs [11]Dp 2.18e-5 cm2/s Diusivity of free O2 in plasma [11]HD 0.45 Hematocrit [11]Jout 0 mmHg cm/s Arteriole outlet 02 diusive ux /α 0 cm/s 02 permeability between arteriole and capillary /

Arteriole vessel wallP0 1.5e-2 cm/s O2 permeability across the lumen-wall interface /Pw 10 mmHg O2 partial pressure in the vascular wall /αw 2.4e-5 mlO2/ml/mmHg [11]Endothelial vessel wallDE 2.8 cm2/s Diusivity of O2 [7]Km,E 150/1.34 mmHg/s O2 consumption rate [7]K1/2,E 4.7 mmHg Michaelis-Menten constant [7]RE −RL 1 µm Thickness of endothelial wall [7]Smooth muscle regionDSMC 2.8 cm2/s Diusivity of O2 [7]Km,SMC 1/1.34 mmHg/s O2 consumption rate [7]K1/2,SMC 1 mmHg Michaelis-Menten constant [7]RSMC −RE 6 µm Thickness of Smooth muscle region [7]

Tissue spacePT 10 mmHg O2 partial pressure in nerve bers layer /hW 2.5e-2 cm/s O2 permeability across the tissue-SMC interface /D 1e-5 cm2/s O2 diusion coecient [3]T 0.025 cm Thickness of the tissue /Km,2 90-180 mmHg/s O2 consumption rate at layer 2 (light-dark) [3]Km,6 26 mmHg/s O2 consumption rate at layer 6 (light-dark) [3]K1/2 2 mmHg Michaelis-Menten constant [3]PC 80-250-405 mmHg O2 partial pressure in Choroid [3]P50 26 mmHg Half partial pressure of O2 saturation in Hb [3]K1/2 2 mmHg Michaelis-Menten constant [3]n 2.7 Hill coecient [3]Pblood 80-250-405 mmHg O2 partial pressure in arterial blood [3]Hb 140 g/L The Hb concentration in blood [3]δ 0.0616 mmol/g O2 carrying capacity of Hb [3]bf 0- 0.2 - 0.4 Blood ow rate in inner retina [3]α1 1.5e-3 nM/mmHg O2 solubility in blood [3]

Table 4.1: Parameters used in the numerical simulation.

75

Conclusions

Very few models in literature address the problem of solving the blood uid-dynamical eld coupled with oxygen transport in each district of the retina.At our knowledge, no models at all present a coupled model, which combinestogether a model for the arterioles, one for the capillary and one for the tis-sue.In this work we have tried to ll this gap, presenting a model for each districtof the retina. We have used an existing model for blood ow in the arterialtree and we have coupled it with the capillary plexi and surrounding tissueusing conservation laws discretized with the Finite Element Method.To this purpose, we have proposed a new geometrical model, based on aDiusion Limited Algorithm, to describe the arterial tree, since in literaturewe only found geometrical structures built as a dichotomic binary system,obtaining a more realistic representation. To show the correctness of theproposed model, we have validated it against several experimental measures.We have also presented novel ideas for future research work. The double-continuum method proposed for the capillary bed may be developed andcoupled with the other models in order to obtain a distribution of hemo-dynamic variables and oxygen partial pressure more detailed than the oneobtained in this work. In this work, we don't have discussed the phenomenonof autoregulation, the intrinsic ability of vascular bed to maintain its bloodow relatively constant despite variations of the arterial pressure Autoregu-lation is pretty important for the correct functioning of mechanisms in theretina and if fails, numerous pathologies can arise, eventually leading toblindness. However, we believe that the models we presented in this workcan be used as a basis for a future research on these phenomena.

76

Appendix A

Calculation of Equivalent

Resistance in Capillaries Bed

In this appendix, we describe the calculation of the equivalent resistancethat we use to treat the coupling between arterial tree and capillaries bed.Assuming that the capillary bed is a continuum, we consider this latter tobe formed by a periodic lattice of 600 unit cells (see Fig.A.1).

Figure A.1: Representation of the capillary bed based on the volume averagingapproach.

Assuming the Hagen-Poiseuille model for uid ow, the capillary do-main can be interpreted as an electric circuit, where the intensity of currentI represents blood ow Q, the potential dierence 4V represents the pres-sure drop 4P while the conductance of the circuit represents the hydraulicconductance of the vessel:

77

Q =πr4

8µl4P ←→ I = G4V. (A.1)

So, we have a circuit as in Fig.A.2

Figure A.2: Capillary bed as an electric circuit. All resistances are equal toeach other.

We assume that all the blood ux coming from the arterial tree is "lumped"at the central point of the capillary bed and there are pressures assigned atthe vertices of this domain. Therefore we insert a current generator at thecentre of the electric circuit and we suppose that the corner vertices of thecircuit are connected to ground. The Iinlet value is arbitrarily chosen.Because of the symmetry of the network, we divide this latter in four squaresand study only one of these. To treat the sides of interface with the othersquare, we replace each resistor R along these with a pair of shunted equalresistors which has an equivalent resistance R (see Fig.A.3)

Using Kirchho's law, i.e. at any node (junction) the sum of currentsowing into that node is equal to the sum of currents owing out of that

78

Figure A.3: A quadrant of the electrical circuit.

node, we nd the value of V in each node of the circuit and in particular V1,the potential value at the node where the current generator is applied. Thenwe can dene the equivalent resistance as:

Req =Vf − V1

Iinlet. (A.2)

Assuming µ = 3.2 ∗ 10−2g/cms, we obtain Req = 1 ∗ 1011g/cm4s.

79

Appendix B

The Finite Element Method

In this Appendix we describe the numerical solution of the model, that wepresented with the 1-D Finite Element Method (FEM)[10]. Let Ω = (a, b) bean open set in R and Γ = ∂Ω the boundary of Ω. The unit outward normalvector to Γ is denoted by −→n and −→v is the given ow velocity.

The general problem consists of nding u = u(x) ∀x ∈ Ω, such that:Lu ≡ ∇ · (−D∇u+−→v u) + σu = f in Ω

u = g on ΓD = x = a(−D∇u+−→v u) · −→n − v+

n u = αu+ β on ΓR = x = b

where D : Ω → R, −→v : Ω → R, σ : Ω → R, f : Ω → R, g : ΓD → R,α : ΓR → R and β : ΓR → R are prescribed data. We consider D ∈ L∞(Ω)with D(x) > 0 a.e., −→v ∈ L∞(Ω), σ ∈ L2(Ω) with σ(x) > 0 a.e., f ∈ L2(Ω),g ∈ H1(Ω), α ∈ L2(ΓR) with α(x) > 0 a.e. and β ∈ L2(ΓR) with β(x) > 0a.e.

We note that the boundary condition on ΓR is equivalent to assigningthe total ux if vn = v−n or to assigning the diusive ux if vn = v+

n .We dene the following function space:

S = u ∈ H1(Ω)| u = g on ΓDV = w ∈ H1(Ω)| w = 0 on ΓD

The objective is to nd u ∈ S, such that

B(w, u) = L(w) ∀ w ∈ V, (B.2)

where

B(w, u) = (∇w,D∇u−−→v u)L2(Ω) + (w, σu)L2(Ω) + (w, u(α+ v+n ))L2(ΓR)

L(w) = (w, f)L2(Ω) + (w, β)L2(ΓR).

Considering a function φ ∈ V and the boundary conditions, the problem(B.2) is equivalent to:

80

Find u ∈ S, such that

φ(b)[α(b) + v+n (b)]u(b)−

∫ b

aJ∂φ

∂xdx+

∫ b

aσuφdx =

=

∫ b

afφdx+ φ(b)β(b) ∀ φ ∈ V.

L : H1(Ω)→ R is a linear and continuous functional, in fact:

|L(w)| ≤‖ f ‖L2(Ω)‖ w ‖L2(Ω) +β(b)w(b)

≤‖ f ‖H1(Ω)‖ w ‖H1(Ω) +C(‖ β ‖L2(Ω)) ‖ w ‖H1(Ω)

≤ (‖ f ‖H1(Ω) +C) ‖ w ‖H1(Ω)

B : H1(Ω)×H1(Ω)→ R is a bilinear and continuous functional:

|B(w, u)| ≤ ‖ ∇w ‖L2(Ω) (‖ D ‖L∞(Ω)‖ ∇u ‖L2(Ω) + ‖ v ‖L∞(Ω)‖ u ‖L2(Ω))+

+ ‖ σ ‖L2(Ω)‖ w ‖L2(Ω)‖ u ‖L2(Ω) +(α(b) + v+n (b))u(b)w(b)

≤ ‖ w ‖H1(Ω)‖ u ‖H1(Ω) (‖ D ‖L∞(Ω) + ‖ v ‖L∞(Ω) + ‖ σ ‖L2(Ω))

+ C(‖ α ‖L2(Ω), ‖ v ‖L∞(Ω)) ‖ u ‖H1(Ω)‖ w ‖H1(Ω)

≤C ‖ w ‖H1(Ω)‖ u ‖H1(Ω)

where C = C(‖ D ‖L∞(Ω), ‖ v ‖L∞(Ω), ‖ σ ‖L2(Ω), ‖ α ‖L2(Ω)),

and coercive functional:

B(w,w) =(D∇w,∇w)L2(Ω) + (σw,w)L2(Ω) + [(α+ v+n )w2]|ΓD

− (vw,∇w)L2(Ω)

but

(vw,∇w)L2(Ω) =

∫Ω∇(vw2)−

∫Ω∇(v)w2 −

∫Ωv∇(w)w

= (vnw2)|ΓR

−∫

Ω∇(v)w2 −

∫Ωv∇(w)w

=1

2[(vnw

2)|ΓR−∫

Ω∇(v)w2]

so that

81

B(w,w) =(D∇w,∇w)L2(Ω) + (σw,w)L2(Ω) + [(α+ v+n )w2]|ΓR

− 1

2(vnw

2)|ΓR+

+1

2

∫Ω∇(v)w2

=(D∇w,∇w)L2(Ω) + (σw,w)L2(Ω) + (αw2)|ΓR+

1

2(∇vw,w)L2(Ω)+

+1

2(vnw

2)|ΓR+

1

2(|vn|w2)|ΓR

− 1

2(vnw

2)|ΓR

=(D∇w,∇w)L2(Ω) + (σw,w)L2(Ω) + (αw2)|ΓR+

1

2(∇vw,w)L2(Ω)+

+1

2(|vn|w2)|ΓR

Hence, the hypotheses of the Lax-Milgram theorem are satised and problem(B.2) has a unique solution.Now we consider a partition of Ω into nite elements ejNj=1 and place

xjN+1j=1 nodes for this partition. Let Sh ⊂ S, Vh ⊂ V be continuous nite

elements spaces where:

Sh = sh ∈ C0(Ω) | sh|ej ∈ P1, ∀j = 1, .., N, sh|ΓD= g;

Vh = vh ∈ C0(Ω) | vh|ej ∈ P1, ∀j = 1, .., N, vh|ΓD= 0.

The classical continuous Galerkin method is:Find uh ∈ Sh, such that

B(wh, uh) = L(wh) ∀ wh ∈ Vh. (B.3)

If we consider φjN+1j=1 as a basis of the space Vh with φj(xi) = δi,j and

ψjN+1j=1 as a basis of the space Vh, we have

wh =

N∑j=1

wjφj and uh =

N∑j=1

ujψj . (B.4)

We set Ωj = x ∈ Ω|φj(x) 6= 0 and we obtain:

j+1∑i=j−1

∫Ωj

(−∂φj∂x

J|Ωj+ σuiψiφj)dx =

∫Ωj

fφj whit j 6= 1, N + 1

82

and

(α(b) + v+n (b))uN+1 −

∫eN

JN∂φN+1

∂xdx+

N+1∑i=N

∫eN

σuiψiφN+1dx =

=

∫eN

fφN+1dx+ β(b).

Using the mass-lumping technique and the Scharfetter-Gummel method (seeAppendix C), the corresponding linear system is:

1 0 0 · · · · · · 0A2,1 A2,2 A2,3 0 · · · 0

0 A3,2 A3,3 A3,4. . . 0

... 0. . .

. . .. . . 0

.... . .

. . .. . .

. . . 00 0 0 0 AN+1,N AN+1,N+1

u1

u2.........

uN+1

=

b1b2.........

bN+1

where

Ai,i−1 = −Di−1

hi−1Be(−vi−1hi−1

Di−1),

Ai,i =Di−1

hi−1Be(

vi−1hi−1

Di−1) +

Di

hiBe(−vihiDi

) + σihi−1 + hi

2, i = 2, ..., N

Ai,i+1 = −Di

hiBe(

vihiDi

),

bi = fihi−1 + hi

2,

and

AN+1,N = −DN

hNBe(−vNhNDN

),

AN+1,N+1 = +DN

hNBe(

vNhNDN

) + α(b) + v+n (b) + σN

hN2,

b1 = g(a),

bN+1 = fN+1hN2

+ β(b).

83

Figure B.1: The red line represents the shape function φt(s).

The nite element method can be used to solve the transport equationin a domain more complex as the one presented to describe the arterial tree.In this case the shape functions have a structure which reects the structureof the domain (see Fig. B.1), i.e. if we consider a node xt and dene withej , ei, ek the only three elements that have xt as node, φt(s) is dened suchthat:

φr(xt) = δt,r,

φt|er (s) = 0 if r 6= j, i, k

φt|er (s) ∈ P1(s) if r = j, i, k.

where s is the curvilinear coordinate.

84

Appendix C

The Scharfetter-Gummel

method

We consider a diusion-transport model in a 1-D domain, Ω = (0, 1), with ahomogeneous Dirichlet boundary conditions:

∂J(u)

∂x= 1 in Ω (C.1a)

J(u) = −µ∂u∂x

+ au with a>0 (C.1b)

u(0) = u(1) = 0. (C.1c)

Dening the weak formulation as:Find u ∈ V , such that:

B(u, v) = F (v) ∀v ∈ V, (C.2)

where

B(u, v) =

∫ 1

0−J(u)

∂v

∂xdx,

F (v) =

∫ 1

01 · vdx,

it can be shown that the hypotheses of the Lax-Milgram theorem are satisedand the problem has a unique solution. Considering a partition Γ = ejNj=1

of Ω, the Galerkin problem associated with (C.2) is :Find uh ∈ Vh(Γ) ⊂ V , such that:

B(uh, vh) = F (vh) ∀vh ∈ Vh, (C.3)

The Peclet number is a dimensionless number relevant in the study oftransport phenomena in uid ows and can be seen as a measure of the

85

relative importance of advection with respect to diusion. Locally, i.e. ineach element ej , it is dened as:

Peloc,j =| a | hj

2µ(C.4)

where hj is the length of element ej . If the Peclet number is larger thatone, the Galerkin solution is corrupted by non-physical oscillation. Whenthe advection term dominates the diusion term in the transport equation,Peloc,j and the Galerkin method loses its best approximation propertyand consequently spurious oscillation appears. To solve this problem, thereare two ways. One is to take h = max | hj | with j = 1, ..., N sucientlysmall, but this can cause a very high computational cost, the other is touse a method that modies the formulation of (C.5) to reduce the Pecletnumber. One of these method is the Scharfetter-Gummel method, in which(C.5) becomes:

Find uh ∈ Vh ⊂ V , such that:

Bh(uh, vh) = Fh(vh) ∀vh ∈ Vh, (C.5)

where

Bh(uh, vh) = B(uh, vh) + bh(uh, vh),

Fh(vh) = F (vh) + fh(vh)

are modied forms with bh and fh proper stabilization terms.We dene

bh(uh, vh) =

∫ 1

0µφ(Peloc)

∂uh∂x

∂vh∂x

dx,

fh(vh) = 0

where φ(t) is a stabilization function, such that:φ(t) > 0 t > 0 (C.6a)

limt→0+φ(t) = 0 (C.6b)

Then we have a new diusion coecient µh,

µh = µ+ µφ(Peloc) ≥ µ

and a new Peclet number Peloc,

Peloc =| a | h2µh

=| a | h

2µ(1 + φ(Peloc))=

Peloc1 + φ(Peloc)

86

that is less than (C.4) i

φ(Peloc) > Peloc − 1. (C.7)

The Scharfetter-Gummel method denes φ(t) as

φ(Peloc) = Peloc − 1 +Be(2Peloc)

whereBe(x) =

x

ex − 1

is the inverse of the Bernoulli function. One can note that if Peloc >> 1,φ(Peloc) ' Peloc and Peloc is less than one.

In particular, if µ and a are constant in each element ej , Jj is constantand is dened as:

Jj = −µjhj

(Be(ajhjµj

)uk −Be(−ajhjµj

)ui) (C.8)

where uk and ui are the values of uh at the nodes of element ej .

87

Bibliography

[1] Bessel Functions of the First and Second Kind. http://www.uma.ac.

ir/files/site1/a_akbari_994c8e8/bessel.pdf.

[2] B. Anand-Aptea and J. G. Hollyelda. Developmental Anatomy of theRetinal and Choroidal Vasculature. 2010.

[3] R. Avtar and D. Tandon. Mathematical Modelling of Intraretinal OxygenPartial Pressure . Trop J Pharm Res, 7(4):11071116, 2008.

[4] S. J. Cringle and D. Y. Yu. A multi-layer model of retinal oxygen supplyand consumption helps explain the mutated rise in inner retinal PO2

during systemic hyperoxia. Comparative Biochemistry and Physiology,132:6166, 2002.

[5] K. Erbertseder, J. Reichold, B. Flemisch, P. Jenny, and R. Helmig. ACoupled Discrete/Continuum Model for Describing Cancer-TherapeuticTransport in the Lung. PLoS ONE, 7(3): e31966., 2012.

[6] P. Ganesan, S. He, and H. Xu. Development of an image Image-Basednetwork model of retinal vasculature . Ann. of Bio. Eng., 38:15661585,2010.

[7] S. I. Gundersen, G. Chen, and A. F. Palmer. Mathematical model ofNO and O2 transport in an arteriole facilitated by hemoglobi based O2

carriers. 143:117, 2009.

[8] A. Harris, G. Guidoboni, J.C. Arciero, A. Amireskandari, L. A. Tobe,and B. A. Siesky. Ocular hemodynamics and glaucoma: the role ofmathematical modeling. Eur. J. Ophthalmol., 23(2):139146, 2013.

[9] R. H. Haynes. Physical basis of the dependence of blood viscosity on thetube radius . 1960.

[10] Thomas J. R. Hughes, Gerald Engel, Luca Mazzei, and Mats G. Larson.The continuous Galerkin Method is Locally Conservative. 163:467488,2000.

88

[11] D. Liu, N. B. Wood, N. Witt, A. D. Hughes, S. A. Thom, and X. Y.Xu. Computational analysis of Transport in the retinal arterial network.Curr. Eye Res., 34:945956, 2009.

[12] B. R. Masters and D. F. Platt. Development of human retinal vessels:embryological and clinical implications. Invest. Ophthalmol. Vis. Sci.,77(30):391, 1989.

[13] D. Moore, A. Harris, D. WuDunn, N. Kheradiya, and B. Siesky. Dys-functional regulation of ocular blood ow: a risk factor for glaucoma?Clin. Ophthalmol., 2(4):849861, 2008.

[14] J. Reichold, M. Stampanoni, A. L. Keller, A. Buck, P. Jenny, and et al.Vascular graph model to simulate the cerebral blood ow in realistic vas-cular networks. Journal of Cerebral Blood Flow & Metabolism, 29:14291443, 2009.

[15] C. E. Riva, J. E. Grunwald, S. H. Sinclair, and B. L. Perring. Bloodvelocity and volumetric ow rate in human retinal vessel. Invest. Oph-thalmol. Vis. Sci., 26:11241132, 1985.

[16] W. M. Roos. Theoretical estimation of retinal oxygenation during retinalartery occlusion. Physiol. Meas., 25:15231532, 2004.

[17] T. F. Sherman. On concerning large vessles to small: the meaning onMurray's law. J. Gen. Physiol, 78:431453, 1981.

[18] C. Simon and I. Goldstein. A new scientic method of identication.New York State J. Medicine, 35(18):901906, 1935.

[19] N. Suwa and T. Takahashi. Morphological and morphometrical analysisof circulation in hypertension and ischemic kidney.

[20] T. Takahashi, T. Nagaoka, H. Panagida, T. Saitoh, A. Kamiya, T. Hein,L. Kuo, and A. Yoshida. A mathematical model for the distribution ofhemodynamic paramters in the human retinal microvascular network.J. Biorheol, 23(7786):29993013, 2009.

[21] P. Tower. The fundus oculi in monozygotic twins:Report of six pairs ofidentical twins. Arch. Ophthalmol., 54:225239, 1955.

[22] F. M. White. Viscous uid ow. 1974.

[23] T. A. Witten and L. M. Sander. Diusion-limited aggregation, a kineticphenomena. Phy. Rev. Lett., 47:14001403, 1981.

89