1 Introduction
We consider two models motivated by the experimental technique known as “hydrogen clearance”, which is used to estimate blood flow in the retina [Reference Alder, Yu, Cringle and Su1, Reference Kety6]. Farrow et al. [Reference Farrow, Hocking, Cringle and Yu4] considered the situation in which the retinal circulation could be considered as a two-layer region consisting of the choroid, in which there is blood flow, and the avascular retina, in which there is not. Here we consider the addition of a third layer in which there is flow on the other side of the avascular retina, and also analyse the effects of variable layer thickness and loss through the exterior of the retina.
The blood in some mammals is supplied to the retina by a combination of the choroidal vascular bed, which is located behind the retina, and the retinal vasculature, which is located within the inner retina. However, some animals, such as guinea pigs, do not have the vascular retina. This extra layer is the difference between the eyes of guinea pigs, who have no retinal vasculature, and humans and rats, who do have a third, vascular layer [Reference Yu, Cringle, Alder, Su and Yu18].
The eye consists of many layers that work to provide clear vision, and it is often considered as a part of the brain [Reference Campbell and Humphries2]. The three major layers in the eye are the outer layer, formed by the sclera (the white of the eye), the cornea (the middle layer, composed of the anterior chamber and the posterior chamber) and the retina [Reference Winchell16]. The choroid is thickest at the extreme rear of the eye at approximately
$ 0.2$
mm, while in the outlying areas, it narrows to approximately
$0.1$
mm [Reference Margolis and Spaide9]. The choroid provides oxygen and nourishment to the outer layers of the retina [Reference Nickla and Wallman12]. The highest oxygen concentration is found in the choroid, and the lowest at the level of the photoreceptor inner segments, reflecting the high oxygen consumption rate in this region. The retina consumes a higher amount of oxygen than any other tissue [Reference Goldman5, Reference Merkle, Leahy and Srinivasan10]. The exchange of metabolites between the blood and the retina is vital to the health of the eye. White et al. [Reference White, Pierce, Nassif, Cense, Hyle-Park, Tearney, Bouma, Chen and de Boer15] found that knowledge of retinal blood flow is an important diagnostic tool in the detection of diseases that lead to blindness such as glaucoma and diabetic retinopathy.
The fundamental vascular bed to supply the nutrition to the retinal pigment epithelium and the outer retina is the choroid, and it plays a vital role in the diseases that infect the eye. Recent studies show that there is a relation between the variations in choroidal thickness and the diseases that affect the eye. For example, with advancing age and myopia, the thickness of the choroid decreases, leading to serious diseases like early age-related macular degeneration [Reference Margolis and Spaide9, Reference Sohn, Khanna, Tucker, Abràmoff, Stone and Mullins13].
There are several techniques to quantify the blood flow through the eye, one of which is the hydrogen clearance technique, which was found to give highly reproducible concentration traces [Reference Alder, Yu, Cringle and Su1, Reference Cringle, Yu, Alder and Su3]. The idea is to inject a bolus of hydrogen saturated saline solution into the blood supply and then measure the hydrogen as it is cleared from the tissue. The principles of hydrogen washout to estimate blood flow are described by Kety [Reference Kety6]. Under ideal conditions, the clearance rate is directly proportional to the local blood flow per unit weight of tissue. The “ideal” conditions require that the measurements should be made in a uniformly perfused three-dimensional tissue, but this situation cannot be guaranteed in the eye due to the short duration of the hydrogen saturated saline injection. Therefore, a model of the situation provides important information to assist in interpretation of the experimental results.
The intraocular blood flow in a rat eye was measured by Yu et al. [Reference Yu, Alder and Cringle17]. Their study used the hydrogen clearance technique to determine retinal blood flow. The arrival and clearance of hydrogen saturated saline to retinal arteries, veins and intervascular were monitored by hydrogen sensitive microelectrodes. They were able to map the time of hydrogen arrival at the retina and removal from the choroid. A mathematical model for the diffusion of inert gases into the retinal tissue and its clearance by blood movement in rabbits was introduced by Strang et al. [Reference Strang, Horton and Gillespie14]. Their study made comparisons between theoretical and experimental data and found that a measure of choroidal blood flow was given by the initial slope of the hydrogen level as it decays. It was also found that the diffusion of inert gases in other ocular tissue can influence the shape of the curve. Farrow et al. [Reference Farrow, Hocking, Cringle and Yu4] presented a simple mathematical model of retinal blood flow that will help to facilitate the study of oxygen supply from the choroid to the avascular retina when there is no vascular retina.
The purpose of this paper is to extend the work of Farrow et al. [Reference Farrow, Hocking, Cringle and Yu4] to consider the effect of the vascular retina. The inclusion of a third vascular layer of variable thickness with leakage of hydrogen through the outer layers will be considered. Some preliminary work was presented by Mansoor et al. [Reference Mansoor, Hocking and Farrow8]. We start with a one-dimensional model of the three layers in which the two outer layers have flow and the middle one does not, and in which there is exchange between the layers with losses through the side. We then move on to consider a two-dimensional model with three layers in a continuum. Note that it is not the purpose of this paper to simulate blood flow in the eye. A description of such a model can be found in the work of Zouache et al. [Reference Zouache, Eames and Luthert19].
2 The three-layer model
This model is a generalization of Farrow et al. [Reference Farrow, Hocking, Cringle and Yu4] who examined the behaviour in the choroid and avascular retina (such as in the retinal system of a guinea pig) [Reference Yu, Cringle, Alder, Su and Yu18]. The three layers to be considered are the choroid, avascular retina and a third layer, the vascular retina, as well-mixed layers (such as in the retinal system of a rat). The blood flow in the choroid and the vascular retina have constant and uniform speeds
$U_{0}$
and
$U_{1}$
, respectively, and there is no flow in the avascular layer. These assumptions are based on the fact that the regions are a mesh of blood vessels rather than a “channel”, so that the flow will be uniform, although in general, one would expect
$U_{0}>U_{1}$
. In addition, we assume there to be exchange coefficients
$\lambda _{0}$
and
$\lambda _{1}$
between the layers and rate of leakage of hydrogen
$\zeta _{0}$
and
$\zeta _{1}$
from the top and bottom layers. This model is the same as the two-layer model (Farrow et al. [Reference Farrow, Hocking, Cringle and Yu4]) if the third layer is removed. Let us define
$\hat C(\hat x,\hat t)$
,
$\hat R(\hat x,\hat t)$
and
$\hat B(\hat x,\hat t)$
to be the concentration of hydrogen in the choroid, avascular retina and vascular retina, respectively, where
$\hat x$
is the distance along the retina and
$\hat t$
is time. A diagram of the main features of the model is shown in Figure 1.
Sketch of the mathematical model in nondimensional variables with layer boundaries
$H_{c}$
,
$H_{ar}$
and
$H_{vr}$
. In the one-dimensional model, each layer is considered as a well-mixed region.
The model equations are, on
$\hat t>0$
,
$0<\hat x<\hat {L}$
,
$$ \begin{align*} \frac{\partial \hat C}{{\partial \hat{t}}}+U_{0}\frac{\partial \hat C}{\partial \hat{x}}=-\frac{{\lambda_{0}}}{\hat{H_{1}}}(\hat C-\hat R)-\frac{\zeta_{0}}{\hat{H_{1}}} \hat C, \end{align*} $$
$$ \begin{align*} \frac{\partial \hat R}{\partial \hat{t}}=\frac{\lambda_{0}}{\hat{H_{2}}} (\hat C-\hat R)+\frac{\lambda_{1}}{\hat{H_{2}}} (\hat B-\hat R) , \end{align*} $$
$$ \begin{align*} \frac{\partial \hat B}{\partial \hat{t}}+U_{1}\frac{\partial \hat B}{\partial \hat{x}}=-\frac{\lambda_{1}}{\hat{H_{3}}} (\hat B-\hat R) -\frac{\zeta_{1}}{\hat{H_{3}}} \hat B , \end{align*} $$
where
$U_{0}$
and
$U_{1}$
are the flow speeds for the choroid and vascular retina, while
${\lambda _{0}}$
and
${\lambda _{1}}$
are the exchange coefficients between these layers and the avascular layer,
$\hat {H_{1}}=\hat H_{c}, \hat {H_{2}}=\hat H_{ar}-\hat H_{c}, \hat {H_{3}}=\hat H_{vr}-\hat H_{ar}$
are the thickness of the choroid, the avascular retina and vascular retina, respectively, and
$\zeta _{0}$
and
$\zeta _{1}$
are the rate of leakage through the outer edges of the choroid and the vascular retina, respectively. The initial conditions are
$$ \begin{align*} \hat{C}(\hat x,0)=f_{0}(\hat x), \quad \hat{B}(\hat x,0)=0 \quad\mbox{and}\quad \hat{R}(\hat x,0)=0, \end{align*} $$
where
$f_{0}(\hat x)$
is the initial concentration distribution of the hydrogen in the choroid, while the other two layers initially have no hydrogen.
The equations can be nondimensionalized using
$\hat {x}=x \hat {H_{1}}$
, so
$L=\hat {L}/\hat {H_{1}}$
is the nondimensionalized length of the region, with
$$ \begin{align*}H_{1}=1, \quad H_{2}=\frac{\hat{H_{2}}}{\hat{H_{1}}}, \quad H_{3}=\frac{\hat{H_{3}}}{\hat{H_{1}}}, \quad \hat{t}=\frac{t \hat{H_{1}}}{U_{0}}, \end{align*} $$
and
$\hat C$
,
$\hat B$
and
$\hat R$
are nondimensionalized using a particular value from the initial condition
$f_{0}(\hat x)$
. Typical values used in this nondimensionalization are
$\hat H_{1} \approx 10^{-4} \approx 100\, \unicode{x3bc} \text {m}$
,
$\hat t = 1$
second and
$\lambda \approx 10^{-7}$
$\text {m}^{2}\,\text {s}^{-1}$
. This means that the length of the region in the simulations corresponds to approximately
$L=3$
mm, and the velocity in the choroid is
$U_{0}\approx 100 \, \unicode{x3bc} \text {m}\,\text {s}^{-1}$
. Then, in nondimensional form,
$$ \begin{align} \frac{\partial C}{\partial t}+\frac{\partial C}{\partial x}=-\alpha_{0}(C-R)- \zeta_{N0} C, \end{align} $$
$$ \begin{align} \frac{\partial R}{\partial t}=\frac{\alpha_{0}}{H_{2}}(C-R)+\frac{\alpha_{1}}{H_{2}}(B-R) , \end{align} $$
$$ \begin{align} \frac{\partial B}{\partial t}+\gamma\frac{\partial B}{\partial x}=-\frac{\alpha_{1}}{H_{3}}(B-R) - \frac{\zeta_{N1}}{H_{3}} B, \end{align} $$
where
$$ \begin{align*}\gamma=\frac{U_{1}}{U_{0}}, \quad \alpha_{0}=\frac{{\lambda_{0}} \hat{H_{1}}}{U_{0}}, \quad \alpha_{1}=\frac{\lambda_{1} \hat{H_{1}}}{U_{0}} ,\quad \zeta_{N0}=\frac{\zeta_{0} \hat{H_{1}}}{U_{0}} \quad \text{and}\quad \zeta_{N1}=\frac{\zeta_{1} \hat{H_{1}}}{U_{0}}\end{align*} $$
are the nondimensional parameters of the problem. The exchange coefficients
$\alpha _{0}$
and
$\alpha _{1}$
are assumed to be constant, as is the ratio of velocities in the two layers. However, it is easy to convert them to be variable parameters if necessary.
If we assume the above equations to be periodic in the domain, then we can write the dependent functions as a complex Fourier series, that is,
$$ \begin{align*} C(x,t)=\sum_{k=-\infty}^{\infty} C_{k}(t)e^{2 \pi i x k/L}, \quad R(x,t)=\sum_{k=-\infty}^{\infty} R_{k}(t) e^{2 \pi i x k/L},\quad B(x,t)=\sum_{k=-\infty}^{\infty} B_{k}(t)e^{2 \pi i x k/L}. \end{align*} $$
Substituting into equations (2.1), (2.2) and (2.3) yields a system of linear ODEs for the coefficients for
$k=-\infty , \ldots ,\infty ,$
$$ \begin{align*} \frac{d C_{k}}{d t}&=-\left(\alpha_{0}+\frac{2 \pi i k }{L}\right)\,C_{k}+\alpha_{0} R_{k}- \zeta_{N0}\, C_{k} , \\ \frac{d R_{k}}{d t}&=\frac{\alpha_{0}}{H_{2}} C_{k}-\left(\frac{\alpha_{0}+\alpha_{1}}{H_{2}}\right) R_{k}+\frac{\alpha_{1}}{H_{2}} B_{k}, \\ \frac{d B_{k}}{d t}&=\frac{\alpha_{1}}{H_{3}} R_{k}-\left(\frac{\alpha_{1}}{H_{3}}+\frac{2 \pi i k }{L}\right)\,B_{k}-\frac{\zeta_{N1}}{H_{3}}\, B_{k}, \end{align*} $$
where initial values of
$C_{k}(0)$
are obtained from the Fourier decomposition of the initial bolus of hydrogen, and
$B_{k}(0)=R_{k}(0)=0.$
We write these equations in matrix form as
$$ \begin{align} \frac { d \hat C_{k}}{dt} =A_{k} \hat {C_{k}} , \end{align} $$
where
$$ \begin{align*}A_{k}=\begin{bmatrix} -\left(\alpha_{0}+\dfrac{2 \pi i k }{L}\right)-\zeta_{N0} & {\dfrac{\alpha_{0}}{H_{2}}} & {0} \\ \alpha_{0} & {-\dfrac{\alpha_{0}+\alpha_{1}}{H_{2}}} &{\dfrac{\alpha_{1}}{H_{3}} } \\ {0} & {\dfrac{\alpha_{1}}{H_{2}} } & {-\left(\dfrac{\alpha_{1}}{H_{3}}+\dfrac{2 \pi i k }{L}\right)-\dfrac{\zeta_{N1}}{H_{3}} } \\ \end{bmatrix} \quad \mbox {and} \quad \hat {C_{k}}=\begin{bmatrix} \ C_{k}\\ \ R_{k} \\ \ B_{k}\\ \end{bmatrix}.\end{align*} $$
This system in equation (2.4) was programmed into the computer package
$\tt {octave}$
and solutions were obtained at different times using the linear multi-step method in the routine “lsode”. The initial condition of
$f_{0}(x)=e^{-x^{2}}/\sqrt {\pi }$
was chosen to represent the initial bolus of hydrogen.
3 Results of the 1-D model
3.1 Constant layer thickness
$\boldsymbol {(H_{1}=H_{2}=H_{3}=1)}$
, no leakage
$\boldsymbol {(\zeta _{N0}=\zeta _{N1}=0)}$
and constant diffusion
$\boldsymbol {(\alpha _{0}=\alpha _{1})}$
When
$\alpha =\alpha _{0}=\alpha _{1}$
, the exchange rate is constant between the three layers. Figure 2 shows a series of profiles of C, R and B at different times when the velocity in the third layer is
$\gamma =0.5$
, the layer thickness is constant, there are no losses through the sides and there is a moderate value of diffusion
$\alpha =0.2$
. The original bolus of hydrogen,
$C(x,0)$
, travels along in the choroid at the speed of flow, slowly diffusing into the other two layers. After a short time, the concentration at the tail end of the original bolus in the choroid drops below the concentration in the avascular retina, resulting in a back diffusion into the choroid. As this occurs, the tail of the original bolus becomes thicker.
This effect can also be seen in the avascular region as there is a return diffusion from the vascular retina. Therefore, as the original bolus moves along, the difference in transfer speed due to the advection and diffusion results in a movement of hydrogen into the avascular retina and then back at the tail end. Figure 2 shows that the original peak value in the choroid is still clear at
$t=12$
, but the maximum value at the tail due to this back flow is becoming larger and at approximately
$t=20$
, it becomes bigger than the leading peak. A similar effect can be seen in the other two layers, as was found by Farrow et al. [Reference Farrow, Hocking, Cringle and Yu4] and Mansoor et al. [Reference Mansoor, Hocking and Farrow8].
This behaviour can clearly be seen in Figure 3 which shows the location,
$x_{m}$
, of the highest value of hydrogen in each of the layers as a function of time for the same parameters as in Figure 2 but for different
$\gamma =0, 0.5, 1$
. The peak value in the choroid seems to travel steadily at the speed of flow until it suddenly jumps backward due to the growth of the secondary peak. The secondary peak (new maximum) in the choroid travels with a new slower speed. Likewise, the peak in the avascular retina starts to travel at the same speed as in the choroid but lags slightly behind, until the peak concentration in the vascular retina leads to a new peak due to the diffusion back from the vascular retina (see Figure 3(a),(b)).
Profiles of
$C(x,t)$
,
$R(x,t)$
and
$B(x,t)$
for
$\alpha =0.2$
at times
$t=4,8,12,16,20$
with
$\gamma =0.5$
. The profiles are shown from top to bottom in order, the choroid (blue), the avascular retina (black) and the vascular retina (red). The avascular and vascular retina profiles have been inverted and the vascular retina has been offset by
$-0.1$
to make the images clearer. (Colour available online.)
Location of the peak concentration,
$x_{m}$
, in the three layers for
$\alpha =0.2$
as a function of time. There is unit velocity in the first layer and
$\gamma =0$
,
$0.5$
and
$1$
in panels (a)–(c), respectively. The choroid is the dashed (blue) line, the avascular retina is the solid (green) line and the vascular retina is the dash–dotted (red) line. (Colour available online.)
In each case, once the transition in the location of
$x_{m}$
has occurred, the new peak value moves at a slower speed. When
$\gamma =0.5$
, the hydrogen in the vascular retina is carried along at half of the speed of the blood velocity in the choroid, which causes the whole process of transition to occur sooner, and also causes the second layer to lag behind the third layer. For
$\gamma =1$
, the process occurs even earlier and very much faster due to the greater mass flux in the third layer. The position of the transition from the current speed of flow to a new speed depends on the values of
$\alpha $
and
$\gamma $
.
Figure 4 shows the location,
$x_{m}$
, of the peak concentration in the three layers for a larger exchange coefficient,
$\alpha =0.5$
, at different values of
$\gamma $
with no leakage and constant layer thickness. This shows the impact of a stronger interchange between the layers. The results obtained for this case show the same qualitative behaviour as for a lower value of
$\alpha $
(see Figure 3), but the jump backwards in peak concentration occurs much earlier due to the more rapid rate of hydrogen flux into the vascular retina. We can clearly see that the shift from the primary peak to the secondary peak in the choroid occurs at
$t\approx 6$
and in the avascular retina at
$t\approx 4$
when
$\gamma =0.5$
. When
$\gamma =1$
, it happens almost immediately in both layers and the first and third layers match in shape because they flow at the same speed, while the second layer lags behind as the motion is dominated by the exchange with the other two layers.
Location,
$x_{m}$
, of the peak concentration in the three layers for a large value of
$\alpha =0.5$
. The velocity of the first layer is
$U_{0}=1$
and in the third layer, it is
$\gamma =0$
,
$\gamma =0.5$
and
$\gamma =1$
(a–c), respectively. The choroid is the dashed line (blue) line, the avascular retina is the solid (green) line and the vascular retina is the dash–dotted (red) line. (Colour available online.)
3.2 Variable layer thickness and diffusion
In this section, we consider the effects of variation in the thickness of the three layers. Typically, the choroid’s thickness is approximately
$\hat H_{1}=100\,{\unicode{x3bc} }\text {m}$
, the avascular retina’s thickness is approximately
${\hat H_{2}=150\,{\unicode{x3bc} }\text {m}}$
and the vascular retina’s thickness is approximately
$\hat H_{3}=50\,{\unicode{x3bc} }\text {m}$
. Here we consider the case with nondimensional widths
$H_{1}=1$
,
$H_{2}=1.5$
and
$H_{3}=0.5$
, since these represent typical values of the ratio of layer thicknesses, as described in [Reference Yu, Cringle, Alder, Su and Yu18]. The results of this modified situation are shown in Figures 5 and 6. Figure 5 shows the location of the maximum hydrogen value in each layer as a function of time for the cases in which
$\alpha _{0}$
is varied while
$\alpha _{1}=0.2$
is held constant and
$\gamma =0,1$
. The figure looks similar to cases with equal layer thickness, but closer inspection reveals that the transition of the maximum point occurs sooner. When
$\alpha _{0}$
takes larger values such as
$\alpha _{0}=1$
or larger, the transition to the terminal velocity occurs almost as soon as the bolus is injected. Therefore, the thicker middle layer (avascular retina) and thinner vascular retina lead to an accumulation of more hydrogen in the middle avascular retina and an earlier peak formation, leading to diffusion back into the choroid, and hence an earlier transition of the maximum hydrogen value. Figure 6 shows the peak location if
$\alpha _{0}=0.2$
and
$\alpha _{1}$
is varied, and
$\gamma =0,1$
. As in Figure 5, the results appear very similar to the cases with equal layer thickness, but the transition in peak concentration occurs earlier but over a slightly shorter time interval.
Location of the hydrogen maximum in the choroid when the layer thicknesses are varied (
$H=1, 1.5\ \mbox {and} \ 0.5$
), respectively, and
$\alpha _{1}=0.2$
is fixed and
$\alpha _{0}$
is varied. In all cases, the varied factor has 6 values which are
$\alpha _{0}$
=
$0.2$
(blue),
$\alpha _{0}$
=
$0.3$
(brown),
$\alpha _{0}$
=
$0.4$
(yellow),
$\alpha _{0}$
=
$0.5$
(purple),
$\alpha _{0}$
=
$1$
(green) and
$\alpha _{0}$
=
$5$
(cyan), at
$\gamma =0$
(a) and
$\gamma =1$
(b). (Colour available online.)
Location of the hydrogen maximum in the choroid when the layers’ thicknesses are varied (
$H=1, 1.5\ \mbox {and}\ 0.5$
), and
$\alpha _{0}=0.2$
is fixed and
$\alpha _{1}$
is varied. In all cases, the varied factor has 6 values, which are
$\alpha $
=
$0.2$
(blue),
$\alpha $
=
$0.3$
(brown),
$\alpha $
=
$0.4$
(yellow),
$\alpha $
=
$0.5$
(purple),
$\alpha $
=
$1$
(green) and
$\alpha $
=
$5$
(cyan), at
$\gamma =0$
(a),
$\gamma =1$
(b). (Colour available online.)
In considering the three layer model in equations (2.1)–(2.3) with different layer thickness, we can see that the effective exchange coefficients include the layer thickness, so if the layer thickness is slightly higher, then diffusion in and out of that layer is slightly reduced (in practical terms, because the travel time across the layer is slightly increased). Therefore, while the initial flow of hydrogen out of the choroid will occur at approximately the same rate, the increase in concentration in the avascular retina will be slightly slower. However, once the level in the avascular retina has peaked, the flow back will occur over a slightly longer time frame, while the value in the choroid will increase quicker. A similar effect will be seen between layers two (avascular retina) and three (vascular retina).
The “terminal” velocity of the peak in the different cases is very similar for each, so that if
$\gamma $
is the same and the layer thicknesses are the same, the terminal values will be the same. This can be seen in the final slopes of the lines in Figures 5 and 6.
Therefore, the effect of different layer thicknesses is similar to changing
$\alpha _{0}$
and
$\alpha _{1}$
in the constant layer thickness case. The effect on the transition is determined largely by the effective value of
$\alpha _{0}$
, while the impact of
$\alpha _{1}$
is to modify the original transition time for the case with
$\alpha _{1}=\alpha _{0}$
by a small amount.
Thus the effective value of
$\alpha _{0}$
(taking into account the variable layer thickness) is the major factor in determining the time of transition to the terminal velocity because it is the controlling factor in the movement of hydrogen from (and back to) the choroid.
Final velocity for different leakage values with two different cases for equal (red dashed line) and variable (blue solid line) thickness, and different velocity in the third layer (
$\gamma =0$
,
$0.5$
and
$1$
), all with
$\alpha _{0}=\alpha _{1}=0.5$
. (Colour available online.)
3.3 Loss through the sides
$\boldsymbol {(\zeta _{N0}, \zeta _{N1} \neq 0)}$
In this section, we address the impact of including the loss of hydrogen through the outer edges of the retina. Figure 7 presents the terminal velocity as the leakage term
$\zeta _{N0}, \zeta _{N1}$
is varied. The case of equal layer thickness is shown as a red dashed line and the blue solid line is for layers of variable thickness, with different values for the velocity in the third layer (
$\gamma =0$
,
$0.5$
and
$1$
). The inclusion of some losses through the edges lowers the level of hydrogen concentration throughout. If the losses are small, say
$\zeta _{N0}=\zeta _{N1}=0.02$
, then not much changes in terms of the transition times and terminal velocity, but if we increase it to
$\zeta _{N0}=\zeta _{N1}=0.1$
, the transition happens earlier than before, while if the leakage term is half of the value of the exchange rate, then the transition will be very early. Blue solid lines present the final velocity with the same parameters but different layer thickness (
$H_{1}=1$
,
$H_{2}=1.5$
and
$H_{3}=0.5$
). It is clear that loss of hydrogen through the sides decreases the level of hydrogen throughout for all situations, and has an impact not only on the diffusion and timing of the transition to the terminal velocity, but also the value of the terminal velocity. Figure 7 also shows a clear impact of variable layer thickness on the terminal velocity as the losses through the side increase, with the variable layer case reaching a lower terminal velocity than the equivalent constant thickness case.
4 Large time behaviour
The results in Section 3 indicate a relationship between the values of velocity in the two outer layers and the final velocity of the peak concentration. An approximate formula for the final velocity of the peak value and also the decay of the hydrogen levels in this model can be derived, and this should assist in untangling the very complicated interactions in this process.
If the system reaches equilibrium as the hydrogen travels along, we equate the terms of the right-hand side from equations (2.1) and (2.3) to give
$$ \begin{align} C=\frac{{\alpha_{0}}\, R} {{\alpha_{0}+\zeta_{N0}}} \quad \mbox {and} \quad B=\frac{{\alpha_{1}}\, R} {{\alpha_{1}+\zeta_{N1}}}. \end{align} $$
Following the analysis of Mitchell et al. [Reference Mitchell, Morton and Spence11], but extending to three layers, adding the three equations (2.1), (2.2) and (2.3) together gives
$$ \begin{align} ({H_{c}} C_{t}+ {H_{ar}} R_{t}+ {H_{vr}} B_{t})+({H_{c}} C_{x}+\gamma {H_{vr}} B_{x})=-{\zeta_{N0}\, C} -{\zeta_{N1} B}. \end{align} $$
By applying the value of C and B from equation (4.1) to equation (4.2), then simplifying the above equation,
$$ \begin{align} R_{t}+V_{\text{term}} R_{x}=-\Gamma_{\text{term}} R, \end{align} $$
where
$V_{\text {term}}$
is the advection velocity of the hydrogen peak and
$\Gamma _{\text {term}}$
is the decay rate:
$$ \begin{align} V_{\text{term}}=\bigg[H_{c} \frac{{\alpha_{0}}} {{\alpha_{0}+\zeta_{N0}}}+\gamma H_{vr} \frac{{\alpha_{1}}} {{\alpha_{1}+\zeta_{N1}}}\bigg]\;\bigg/\;\bigg[H_{ar}+ \frac{H_{c} {\alpha_{0}}} {{\alpha_{0}+\zeta_{N0}}}+ \frac{H_{vr} {\alpha_{1}}} {{\alpha_{1}+\zeta_{N1}}}\bigg] \end{align} $$
and
$$ \begin{align} \Gamma_{\text{term}}=\bigg[ \frac{ \zeta_{N0} {\alpha_{0}}} {{\alpha_{0}+\zeta_{N0}}}+ \frac{ \zeta_{N0}{\alpha_{1}}} {{\alpha_{1}+\zeta_{N1}}}\bigg]\;\bigg/ \;\bigg[H_{ar}+ \frac{H_{c} {\alpha_{0}}} {{\alpha_{0}+\zeta_{N0}}}+ \frac{H_{vr} {\alpha_{1}}} {{\alpha_{1}+\zeta_{N1}}}\bigg]. \end{align} $$
Some results obtained from equations (4.4) and (4.5), applied to different cases, are given in Table 1. These values agree very well with those in the three-layer simulations shown above. Figure 7 confirms the general behaviour of the final velocity values for different velocities
$\gamma $
in the vascular retina with
$\alpha =0.5$
with constant and different layer thickness, and different rates of leakage,
$\zeta $
. In all cases, all the parameters have a different impact on the final velocity, for instance, when the thickness of the layers is constant and there is no leakage, different
$\alpha $
does not effect the final velocity of the peak concentration. However, if we include these effects,
$\alpha $
will have an impact on the final velocity. Figure 8 shows the decay from the simulations compared with the value
$e^{-\Gamma _{\text {term}}t}$
for the peak value of concentration, and shows excellent agreement with the predicted values. Therefore, equations (4.3), (4.4) and (4.5) provide an excellent indicator of all of the behaviours identified in the preceding sections, and the influence of flow speed, layer thickness, losses and diffusion rates in this model.
Final velocity values for different parameters
$\alpha $
,
$\gamma $
,
$thickness\ H$
and
$\zeta $
.
Decay rate in the retina given from
$\Gamma _{\text {term}}$
in equation (4.5) at different values of leakage
$\zeta =0.1$
(a) and
$\zeta =0.3$
(b) with
$\alpha =0.5$
and constant layers thickness H. Regression curves for peak values of choroid,
$\exp (-0.07 t)$
for
$\zeta =0.1$
and
$\exp (-0.19t)$
for
$\zeta =0.3$
. (Colour available online.)
5 Two-dimensional model
To study this problem more deeply, we can extend the model to be two-dimensional and treat it as a single layer including the transport by diffusion across all of the three layers (the choroid, the avascular retina and the vascular retina) (see Figure 1). Therefore,
$$ \begin{align*} \frac{\partial C}{\partial {t}}+U(z)\frac{\partial C}{\partial {x}}=\frac{1}{Pe} \left(\frac{\partial^{2} C}{\partial {x^{2}}}+\frac{\partial^{2} C}{\partial {z^{2}}}\right), \end{align*} $$
where
$C(x,z,t)$
is the hydrogen concentration over all of the three regions,
$Pe=U_{0} \hat {H_{1}}/\kappa $
is the Peclet number,
$\kappa $
is the coefficient of diffusion in the eye tissue and
$U(z)$
is the nondimensional velocity,
$$ \begin{align*} U(z) = \begin{cases} 1 & \text{for } -1< z\leq 0,\\ 0 & \text{for } -H_{AR}< z\leq -1,\\ \gamma & \text{for } -H_{VR}\leq z\leq -H_{AR}, \end{cases} \end{align*} $$
where
$H_{AR}=1+H_{2}$
and
$H_{VR}=1+H_{2}+H_{3}$
are the boundary locations for the choroid, the avascular retina and the vascular retina, respectively, in nondimensional terms. The quantity
$\kappa $
is analogous to and of a similar magnitude to the values of
$\lambda $
in the three-layer model. Here we have set all three layers to potentially be of different thickness.
To consider the losses through the sides of the retina, the boundary conditions are
$$ \begin{align*} \frac{\partial C}{\partial t}=\zeta_{L} C \quad \mbox{ on } z=0,-H_{vr}, \end{align*} $$
where
$\zeta _{L}$
is the leakage rate, so that zero flux corresponds to
$\zeta _{L}=0$
. The initial condition consists of a sine profile bolus of hydrogen in the choroid, i.e.
$$ \begin{align*} C(x,z,0) = \begin{cases} -\sin(\pi x) \sin(\pi z), & \text{for } \;\,0\leq x\leq 1, -1\leq z\leq 0,\\\\ 0,& {\mbox {otherwise}}. \end{cases} \end{align*} $$
The two-dimensional equations were solved numerically using a finite volume method. Quadratic upwinding was used for the advection term [Reference Leonard7] and a standard centred difference was used for the diffusion term. The resulting large system of ordinary differential equations was solved by using routine
$lsode$
from the computer package
$\tt {octave}$
, to step through time. Note that the Peclet number behaves as the inverse of the coefficient of diffusion in the three-layer model. Large Peclet number, therefore, means small diffusion.
5.1 Simulation results for the 2-D model
Many simulations were conducted for different values of
$Pe$
, thickness and leakage. However, we will only discuss the key results. Figures 9, 10 and 11 show contours of the concentration of hydrogen for a large Peclet number,
$Pe=12$
, at different times for different velocity ratios,
$\gamma $
, between the layers in the case with equal layer thickness and no leakage. This value of
$Pe$
corresponds to relatively low diffusion.
Distribution of hydrogen in the three layers for
$t=0,10, 30$
and
$50$
and
$Pe=12$
at
$\gamma =0$
with equal layer thickness and no leakage. Contours at the right-hand end at
$t=30$
are due to the periodic boundary condition. (Colour available online.)
Distribution of hydrogen in the three layers for
$t=30$
and
$50$
and
$Pe=12$
at
$\gamma =0.5$
with equal layer thickness and no leakage. Times
$t=0,10$
look very similar to those in Figure 9. (Colour available online.)
Distribution of hydrogen in the three layers for
$t=30$
and
$50$
and
$Pe=12$
at
$\gamma =1$
with equal layer thickness and no leakage. Times
$t=0,10$
look very similar to those in Figure 9. (Colour available online.)
For the case with zero velocity in the vascular retina, i.e.
$\gamma =0$
and
$t=0$
and
$t=10$
(Figure 9), the bulk of the hydrogen remains in the choroid, while the vertical spread across the region occurs slowly. Thus advection in the choroid “drags” the flow along the horizontal region with a trailing bolus of hydrogen through the other two layers, so that there is a maximal region stretching from a long way along the choroid to a lot further back in the vascular retina (which in this case has no flow). The contours at times
$t=0$
and
$t=10$
are very similar in the cases for Figures 10 and 11, since the hydrogen has not yet reached the vascular retina, and so are not shown.
This trailing effect is reduced in the case of
$\gamma =0.5$
(Figure 10), but is still present as a distributed region of higher concentration extending backward and downward, but not as far. If
$\gamma =1$
(Figure 11), this affect is mitigated even further, and the maximal region curves from the choroid backward into the avascular retina and then forward again into the vascular retina. The difference is clearly due to the differing transport in the third, vascular retina layer, so that as
$\gamma $
increases, the work of carrying the hydrogen along is conducted more by the third layer.
When
$t=30$
and
$t=50$
with
$\gamma =0$
,
$\gamma =0.5$
and
$\gamma =1$
, we notice that the peak hydrogen region is spread across the full width. This is evidence of the equivalent of the slowing of the advection of the peak hydrogen concentration observed in the previous sections. In the case of
$\gamma =1$
(Figure 11), the structure becomes symmetric because the flow velocity in the two outer layers is the same. However, the highest concentration always remains in the choroid, and we can conclude that the advection term dominates, as the large value of the Peclet number is directly comparable to the small value of
$\alpha $
in the three-layer one-dimensional model.
Figure 12 demonstrates the comparable contours to those in the previous figures of the distribution of hydrogen, but with small
$Pe=1$
and
$\gamma =1$
. The diagrams for
$\gamma =0, 0.5$
are not shown, but look very similar. For this value of
$Pe$
, the diffusion term plays a vital role in the transport of the hydrogen, with the spread across the layers being very rapid. As a result, the contours appear quite similar for all cases of
$\gamma =0, 0.5, 1$
with almost vertical contours around the location of the original bolus of hydrogen travelling forward.
Distribution of hydrogen in the three layers for
$t= 10$
and
$30$
and
$Pe=1$
at
$\gamma =1$
with equal layer thickness and no leakage. The hydrogen has spread evenly across all layers. (Colour available online.)
Figure 13 show contours analogous to Figure 11, with
$Pe=12, \gamma =1$
but with different layer thicknesses as considered in the one-dimensional model. We notice that we get similar results to the case where the layers have equal thickness, but with a delay in vertical spread across the region and a loss of symmetry when
$\gamma =1$
. When
$Pe=1$
, the diffusion becomes dominant and the effect of layer thickness is greatly reduced, and so contours for this case look very much like those in Figure 12. The impact of variable layer thickness can be seen clearly in the value of the final velocity, such as in Figure 13 with
$\gamma =1$
. This figure is similar to the results for
$\gamma =0.5$
when the layers are of equal thickness, since the thinner layer results in less advective transport.
Distribution of hydrogen in the three layers of variable thickness for
$t=30$
and 50 and
$Pe=12$
at
$\gamma =1$
with no leakage,
$\zeta _{L}=0$
. Contours at times
$t=0,10$
look very similar to those in Figure 9. (Colour available online.)
When
$Pe=1$
, the speed of travel of the peak hydrogen concentration rapidly transitions to the slower, terminal velocity from equation (4.4). Noting that
$Pe=1$
is equivalent to a large value of
$\alpha $
, this is consistent with the earlier transition in velocity seen in the one-dimensional model.
Figure 14 displays the profile of the hydrogen with a large value of
$Pe=12$
in the case of the velocity in the vascular retina
$\gamma =0.5$
with different layer thickness and a small value of leakage
$\zeta _{L}=1$
. An interesting result for this case is that the location of the bulk of the hydrogen in the choroid has been affected by including the leakage term. At
$t=10$
, the bolus of the hydrogen starts moving along the choroid and spreading to the other layers but the bulk is still in the choroid. However, by
$t=20$
and
$t=30$
, the hydrogen is moving along the choroid very slowly due to the different layer thickness and the leakage. Therefore, the location of the bulk of the hydrogen shifts into the avascular retina.
Distribution of hydrogen in the three layers for
$t=10, 20$
and 30 and
$Pe=12$
at
$\gamma =0.5$
with different layer thickness including small leakage rate
$\zeta _{L}=1$
. (Colour available online.)
Figure 15 displays the profile of the hydrogen with a large value of
$Pe=12$
in a case where the velocity in the vascular retina is
$\gamma =0.5$
with different layer thickness including a large value of the leakage term
$\zeta _{L}=6$
. Similar results to Figure 14 are obtained, but the development of the maximum in the avascular retina happens much quicker, so that it has occurred by
$t=10$
. In addition, the peak concentration in the choroid is travelling at slower speed than that set by the advection velocity. The peak concentration at
$t=20$
is located at
$x\approx 12.5$
and at
$t=30$
is located at
$x\approx 15.5$
, whereas under pure advection, it would be at
$x=20$
and
$x=30$
. This is further evidence of the velocity transition seen in the one-dimensional model.
Distribution of hydrogen in the three layers for
$t=10, 20$
and 30 and
$Pe=12$
at
$\gamma =0.5$
with different layer thickness including large leakage rate
$\zeta _{L}=6$
. (Colour available online.)
Figure 16 shows the concentration of hydrogen with a large value of
$Pe=12$
when the velocity in the vascular retina
$\gamma =1$
and with equal layer thickness and a large value of leakage
$\zeta _{L}=6$
. The results are different to those in Figure 15 because the layers are now of the same thickness. A symmetric profile develops very rapidly and then remains as it travels along.
Distribution of hydrogen in the three layers for
$t=10, 20$
and 30 and
$Pe=12$
at
$\gamma =1$
with equal layer thickness including large leakage rate
$\zeta _{L}=6$
. (Colour available online.)
Figure 17 shows the contours of the hydrogen with a small value of
$Pe=2$
when the velocity in the vascular retina
$\gamma =1$
with different and equal layer thickness, including large values of leakage
$\zeta _{L}=1$
at
$t=30$
. The pictures are very similar except that there is a slower travel speed and a slight asymmetry in the case of variable layer thickness because of the lesser flux in the vascular retina. The faster rate of diffusion for
$Pe=2$
means that the advection dominates the flow once the hydrogen has spread across the whole retina, and the leakage again leads to the centre of the hydrogen being situated in the avascular retina.
Distribution of hydrogen in the three layers for
$t=30$
and
$Pe=2$
at
$\gamma =1$
with different (a) and constant (b) layer thickness including large leakage rate
$\zeta _{L}=1$
. (Colour available online.)
6 Final comments
We have considered two models of the flow of hydrogen in the retina relevant to the method of estimating retinal blood flow known as the “hydrogen clearance method”. The work of Farrow et al. [Reference Farrow, Hocking, Cringle and Yu4] is extended to consider a third “vascular retinal” layer to take account of the differences for example between the eyes of rats and guinea pigs [Reference Yu, Cringle, Alder, Su and Yu18], as well as different layer thickness and leakage through the sides. The work is not meant to be a simulation of blood flow in the retina such as in Zouache et al. [Reference Zouache, Eames and Luthert19], which is a much more complicated problem, but simply to consider the effects of layered advection and diffusion. There are a number of approximations inherent in these models. In reality, the layers do not run for long distances in a single direction and the vascular layers are not single flows but complicated meshes of blood vessels. However, these assumptions allow us to see the interesting interactions between the advection, diffusion and loss terms unhindered by complicated geometry.
The presence of the third, vascular retinal layer has been shown to have an effect on the value of terminal advective velocity, the timing of the transition to this velocity and also on the shape of the concentration contours within the retina. Knowledge of this third layer is therefore important in determining the flow of hydrogen, and hence for the interpretation of the hydrogen clearance curves. Interestingly, as the rate of loss through the sides increases, the centre of concentration of the hydrogen shifts from the choroid into the middle of the avascular retina. The time to transition to the slower, terminal advection velocity is determined much more by the diffusivity than the speed of flow, and occurs earlier for larger rates of diffusion (smaller Peclet number). In particular, the terminal advection velocity given by equation (4.4) and decay rate of hydrogen concentration given by equation (4.5) can be obtained with great accuracy from within the new large time equation (4.3). This equation provides an excellent guide for comparison with experimental values, and is a very good representation of the process as time goes on. Comparison with the experimental results suggests that the Peclet numbers in the retina are more in the lower range of
$Pe=1,2$
than higher values, meaning that the transition to the lower, terminal velocity will occur quite soon after the injection of the bolus of hydrogen, a fact that needs to be considered when examining the experimental results.
The results of the well-mixed, layered model are qualitatively confirmed by the two-dimensional, continuous model, which showed the manner in which the hydrogen spread across the whole region. In this simpler model, it is shown that variable layer thickness can be accommodated by modification of the diffusivity and flux in the vascular retina layer. Together, the two models provide a very good analysis and understanding of the spread of hydrogen within this representation of the retina.
