Aqueous humor dynamics in human eye: A lattice Boltzmann study

1.   Introduction

Glaucoma is a severe progressive ocular disease, which can result in vision loss and is the 2rd most common cause of blindness worldwide ^[1]. As elderly population grows, more and more people are expected to be affected by this disease. For a long time, prophylaxis and treatment of glaucoma have attracted substantial attentions of the researchers all over the world. Lots of researches have shown that glaucoma is closely related with the dynamics of aqueous humor (AH) ^[2]. The anterior segment of eye is filled with AH, a clear, colorless, water-like fluid that is continuously secreted by the ciliary body (CB), located just posterior to the iris. It passes through the gap between the iris and the lens into the anterior chamber (AC), and then circulates there due to natural convection. Eventually, the bulk of AH drains from the eye via specialized tissues consisted of the trabecular meshwork (TM), the juxtacanalicular connective tissue (JCT), the endothelial lining of Schlemm's canal (SC), the collecting channels and the aqueous veins ^[3]. As these specialized tissues can cause a significant hydrodynamic flow resistance, a positive pressure in the eye is required to drain the AH out of the eye, and the positive pressure is called intraocular pressure (IOP). IOP takes a central role in the AH dynamics, which stabilizes the otherwise flaccid eye and ensures its normal physiological functions. When an obstruction of TM or SC occurs, an elevation in IOP is then observed, and this direct mechanical effect is the most common risk factor for glaucoma progression. High-level of IOP, sustained for a long time, can damage the optic nerve in the eye and eventually result in vision loss, which is the primary cause of blindness in glaucoma. Fortunately, several studies have shown evidence that IOP can be controlled by decreasing the production of AH and/or decreasing the hydrodynamic resistance of the trabecular outflow pathways for the drainage of AH ^[4,5,6]. Therefore, revealing the pathogenesis and development mechanisms of glaucoma, in terms of the dynamics of AH, is crucial to its prophylaxis and therapy. Over the years, substantial efforts have been devoted to better understanding the dynamics of AH, both analytically and numerically ^[2,7].

Glaucoma is often, although not always, accompanied by an elevation in IOP, and this elevation is usually known to be resulted from an increase in the hydrodynamic resistance to AH drainage. Originally, lots of work focused on the analysis of TM as a porous material filled with a gel-like biopolymer ^[8], and the predictions on the overall permeability of TM by this method is quite close to the experimental data. However, the subsequent treatments of TM, according to the predictions of this model, showed that the effect of reducing hydraulic resistance is much smaller than would be expected ^[9]. In another hand, it is also important for the mechanism of AH, which drains from the eye via the small pores in the cellular lining of Schlemm's canal ^[10]. Assumed the flow through a single pore and modeled the whole cellular lining of SC as a parallel network consisting of hydrodynamically non-interacting pores, the resistance of a single pore to AH outflow can be estimated by Sampson's theory in low-Reynolds number hydrodynamics. However, the results are much lower than the observed total resistance of aqueous outflow ^[2]. Johnson et al. ^[11] noted a fact that there is a hydrodynamic interaction between the upstream porous material of TM and the endothelial lining of SC. Based on this fact, a simplified unit-cell model, composing of a single pore and the upstream region of the porous medium drained by that pore, was proposed, with which the calculated hydraulic resistance of the ensemble structure is in agreement with experimental measurements ^[12]. By using the finite element method, Bradley and Heys ^[13] researched the effects of variable permeability on AH outflow, and found that the overall resistance of JCT would be increased when the tissue is assumed to be heterogeneous (instead of homogenous).

In addition to the drainage of AH from the TM, the overall flow process of AH in human eye, involving AH generation, flow, circulation, and finally drainage, attracts researchers' attentions all the way. By using analytic method, Avtar and Srivastava ^[14] constructed a simple mathematical model by a series of approximations, discussed the natural convectional flow of AH in the AC of a human eye, and gained the expressions for the temperature and velocity profiles arising from the temperature gradient across the AC. Considering the coupling of Navier-Stokes and Darcy flows, Crowder and Ervin ^[15] proposed a numerical method for investigating the fluid pressure in human eye without the temperature difference across the AC, and found that the increase in IOP is related to the viscosity of AH, the permeability of TM and the flow rate of AH. In order to investigate the drug delivery through the cornea, from a therapeutic lens to the AC of eye, Ferreira et al. ^[6] established a mathematic model for the AH flow in both healthy and pathologic eyes, and numerically studied the interaction of drug flow with the dynamics of AH. The results showed that the drug delivery through the lens has an obvious delay. In consideration of the existence of segmental outflow, Loke et al. ^[16] developed a numerical model for investigating the influences of the segmental AH outflow and the eye orientation on ocular drug delivery, and they suggested the design of ophthalmic drug delivery is needed to re-evaluate according to these effects. Using the method of theoretical analysis combined with numerical calculation, Fitt and Gonzalez ^[17] studied in detail the flows of AH in the AC driven by various physical mechanisms, and they found that the buoyancy caused by the intraocular temperature gradient is the most common mechanism leading to the AC flow, generating velocities are orders of magnitude faster than those generated by other physical mechanism. Zhang et al. ^[18] proposed a 3D computational model for drug delivery in anterior eye after subconjunctival and episcleral implantation, and found that subconjunctival implantation is more effective than episcleral implantation. In addition, Canning et al. ^[19] analyzed the flows of AH in the AC with/without particles by the theoretical and numerical methods, and found that their models can well predict the established and observed features that may exist in a traumatized eye. Heys et.al ^[20] presented a mathematical model of the anterior eye, and explored the elastohydrodynamic effects of accommodation on both the contour of the iris and the AH dynamics. The results confirmed that accommodation produces bowing of the posterior iris and the magnitude of the bowing is a strong function of the amount of accommodation. In order to investigate the heat transfer in the human eye, Karampatzakis and Samaras ^[21] presented a 3D numerical model, and found that the consideration of AH fluid dynamics would affect the temperature distributions on the corneal and lens surfaces.

In conclusion, although great progress has been made in the study of AH dynamics in human eye. However, there still are many unrevealed knowledge and challenges about the AH dynamics to be further explored ^[2]. At the same time, as far as we know, the studies on the AH dynamics are generally for a healthy eye ^{[14,15,16,17,18,19,20,21]}, while few systematic studies on the AH dynamics are for a pathologic eye. The pathological changes of intraocular tissue may have an important effect on the AH dynamics, and some meaningful results may be obtained by studying them. In view of the methods, all above mentioned works have been completed either by theoretical methods or numerical ones. For the theoretical methods, as a series of approximations are necessary, the methods have limited applications in some special circumstances. As a potentially promising method, the numerical simulation has been successfully applied to various fields. However, as the inherent complex anatomical structure of the human eye, together with the high diversity of the tissues, the aqueous humor flow is a typical complex one, and to more clearly understand the AH dynamics is confronted with dire challenges. Therefore, developing an effective model with simple algorithm, easily dealing with complex geometric boundaries and higher accuracies is significant.

Lattice Boltzmann method (LBM), as an alternative new computational hydrodynamics method, has achieved great progress in recent years ^{[22,23,24,25,26,27]}. The inherent advantages of simple algorithm, easy implementation of complex boundaries, high efficiency and full parallelization enable its successful applications on modeling complex fluid systems ^{[28,29,30,31,32]}, including those in industry, science, bioengineering, biomedical engineering and so on.

In view of the dimension of human eye and the complex flow involving complex boundaries and coupling of the multi-physics problems, we can expect that LBM, a newly developed mesoscopic numerical method, may be a good candidate for modeling the dynamics of AH. However, to our best knowledge, there have been few reports of modeling the dynamics of AH with LBM. Thus it is quite natural that, this paper will focus on developing a coupled lattice Boltzmann model for the dynamics of AH in the human eye, and then checking its validity. By using the proposed model, we want to systematically investigate the roles of various macroscopic quantities, such as inlet flow rate, viscosity of AH, permeability of TM, anterior chamber angle and so on, in the dynamics of AH both for healthy and pathologic eyes, as well as their influences on IOP, and expect to obtain some interesting results.

The paper is organized as follows. First, in section 2, we describe the development of a coupled lattice Boltzmann model for the dynamics of AH, including the macroscopic governing equations and the corresponding lattice Bhatnagar-Groose-Krook models. Section 3 starts to check the validity of the proposed model through a case of healthy eye, and then several additional cases are conducted in both healthy and pathologic eyes to explore the effects of various parameters, such as the inlet velocity, viscosity of AH, permeability of TM, anterior chamber angle, etc. on IOP, and after that, the wall shear stresses on the corneal endothelium and the iris surface are investigated. Finally, we draw some conclusions.

2.   Lattice Boltzmann model for the dynamics of AH

The dynamics of AH occurs in the anterior segment of eye, made up of cornea, TM, conjunctiva and sclera externally, and consisting of AC, iris, pupil, posterior chamber (PC) and ciliary body(CB) internally. It is characterized with behaviors of a multi-scale multi-physics coupled system.

2.1.   Geometry

The anterior segment of eye consists of two chambers: the anterior and posterior chambers filled with AH. Based on the anatomy and physiological dimensions of the human eye ^[33], a 2D geometry resulting from a horizontal section is given in Figure 1, the posterior and anterior chambers of the eye is represented by $\Omega_{1}$ . The section of the iris is expressed as $\Omega_{2}$ , and the section of TM, an annular structure (in the basis of the cornea), is denoted as $\Omega_{3}$ . The interface between the cornea and AC is represented by $\Gamma_{1}$ , and the bottom surface between the lens and PC as $\Gamma_{2}$ . All these sections and the interfaces are signaled in Figure 1 and the related geometric dimensions are listed in Table 1.

2.2.   Governing equations

To simulate the dynamics of AH in the anterior segment of eye, we firstly take into account that the AH is a transparent and clear liquid with almost the same viscosity as saline water, can be treated as Newtonian. The cornea, iris and lens are assumed to be solid materials, without any changes in shape during AH flows, while the TM is thought as an isotropic, homogenous, porous matrix.

2.2.1.   Incompressible flow in the anterior and posterior chambers

Once aqueous humor is secreted from the ciliary body, at a flow rate of 2 to 2.5 μL/min, it flows forward to the posterior chamber, through the pupil into the anterior chamber, and then circulates there due to natural convection. In view of the very small flow rate and modest dimension, the flow of AH is creeping, compressibility and inertia can be neglected. Therefore, without considering the heat transfer effect, the flow of the AH in the posterior and anterior chambers can be described with the incompressible Navier-Stokes equations.

where $\boldsymbol{u}$ is the velocity of AH, $\nu$ kinematic viscosity, $p$ the pressure in AH, and $\boldsymbol{F}$ is the body force. In order to solve Eqs (1) and (2), no-slip boundary conditions are imposed on the surfaces ${\mathrm{\Gamma }}_{1}$ , ${\mathrm{\Gamma }}_{2}$ and that of iris, while velocity boundary conditions are applied on the AH inlet (at ciliary body) and outlet (at trabecular meshwork) boundaries.

2.2.2.   Convection and diffusion of heat in anterior chamber

In this investigation, the effects from the convection and diffusion of heat, induced by the temperature difference cross the AC, on the dynamics of AH are in consideration. In practice, there would be a temperature difference between the inner surface of the cornea, ${\mathrm{\Gamma }}_{1}$ , and that of the lens, ${\mathrm{\Gamma }}_{2}$ , thus the convection and diffusion of heat in the AC is inevitable. It is known that if the viscous heat dissipation and compression work done by the pressure can be ignored, the temperature field is passively advected by the fluid flow and obeys a simple passive-scalar equation ^[34].

where $T$ is the temperature in the temperature field, and $\theta$ the diffusivity. The heat flux boundary conditions with specific temperatures are imposed on the inner interfaces of the anterior and posterior chambers, which are coupled to Eq (3).

2.2.3.   Seepage of aqueous humor trough trabecular meshwork

The TM is assumed as a cylindrical annular ring surrounding the AC with specific thickness, approximatively described as an isotropic, homogeneous, porous matrix swollen with continuously flowing AH in this work. According to the porous media theory, the seepage of AH obeys the Darcy's law, which can be written as

where $\boldsymbol{u}$ is the velocity of the seepage, $p$ the pressure in the TM, $\mathrm{K}$ hydrodynamic permeability, and $\mu$ the viscosity of AH. The pressure boundary conditions with specific pressures are applied on the two sides of TM.

2.2.4.   Coupling of multiple physics

In the present dynamics of AH, incompressible Navier-Stokes flows, convection and diffusion of heat and Darcy's flow are coexistence, which should be coupled with each other. As the flow of AH in the anterior and posterior chambers is natural convection, the well-known Boussinesq approximation can be used to couple the incompressible Navier-Stokes flow with the convection and diffusion of heat. In the Boussinesq approximation, the physical parameters of fluid such as density, viscosity and thermal diffusion coefficient are assumed to be constant except the density in external force term, where the fluid density $\rho$ is assumed to be a linear function of the temperature

where ${\rho }_{0}$ and ${T}_{0}$ are the average fluid density and temperature, respectively, $\beta$ is the thermal expansion coefficient.

For the incompressible Navier-Stokes flow and the Darcy's flow, they are slightly simpler to handle. At first the IOP is calculated with Eq (2), and then the velocity of the seepage can be obtained through the difference between the IOP and the pressure of the vein, ultimately this desired velocity will be used as the velocity boundary condition, enabling the iteration of Eq (2) to be repeated continually. With this algorithm, the coupling of Eq (2) with Eq (4) is implemented.

2.3.   Lattice Boltzmann model for dynamics of AH

In lattice Boltzmann method, we use density distribution functions ${f}_{i}\left(\boldsymbol{x}, t\right)$ to depict a fluid, and the distribution functions are represented the probability that a pseudo-fluid particle with velocity ${\boldsymbol{e}}_{i}$ comes into the lattice site $\boldsymbol{x}$ at time $t$ . In order to mimick the motion of the pseudo-particles, at each time step, the pseudo-particles entering the same lattice site collide, and then the resulting distribution functions are streamed to the neighboring sites. The admissible velocities ${\boldsymbol{e}}_{i}(i = \mathrm{0, 1}, \dots z)$ , of component ${e}_{i\alpha }$ , are dependent on the lattice topology, where z represents the lattice coordination number (i.e., the number of lattice links). Conventionally, ${\boldsymbol{e}}_{0} = 0$ and ${f}_{0}\left(\boldsymbol{x}, t\right)$ is the density distribution of particles at rest ^[35]. In single-relaxation-time approximation, the distribution functions obey the following evolvement equation in the lattice unit

with $\delta x = c\delta t$ is the lattice spacing, $c$ the lattice velocity (magnitude), $\delta t$ the time step, and ${F}_{i }$ is a body force term if it exists. In the lattice Boltzmann dynamics of AH, it is obvious that the velocity field and the temperature field should be solved respectively by using two independent lattice Bhatnagar-Gross-Krook (BGK) equations, and then combined into a coupled system with the Boussinesq approximation.

As for the velocity field, the LBM for incompressible flow proposed by Guo et al. ^[34] is adopted, in which the equilibrium distribution function takes the form

where $\sigma$ ¸ $\lambda$ and $\gamma$ are parameters satisfying $\lambda +\gamma = \sigma$ and $\lambda +2\gamma = 1/2$ . In Eq (7), ${S}_{i}\left(\boldsymbol{u}\right)$ is a function of macroscopic velocity $\boldsymbol{u}$ and discrete velocity ${\boldsymbol{e}}_{i }$

with the weight coefficients ${\omega }_{0} = 4/9$ , ${\omega }_{1} = {\omega }_{2} = {\omega }_{3} = {\omega }_{4} = 1/9$ and ${\omega }_{5} = {\omega }_{6} = {\omega }_{7} = {\omega }_{8} = 1/36$ for D2Q9 topology. The primitive macroscopic variables of the incompressible fluid, the velocity $\boldsymbol{u}$ and pressure $p$ , are given by

Through a multi-scaling expansion procedure, the incompressible Navier-Stokes equations (1) and (2) can be recovered from this incompressible lattice BGK model to the order of $O\left(\delta {t}^{2}\right)$ or $O\left(\delta {x}^{2}\right)$ , if the microscopic velocity $c = O\left(1\right)$ , in which the kinematic viscosity $\nu$ is given by

The corresponding no-slip boundary conditions on the surfaces of ${\mathrm{\Gamma }}_{1}$ , ${\mathrm{\Gamma }}_{2}$ and that of iris are implemented with the method proposed by Bouzidi et al. ^[36], and the velocity boundary conditions, on the inlet (ciliary body) and the outlet (trabecular meshwork), are evaluated by the algorithm presented by Guo et al. ^[37].

In regard to the temperature field, a lattice with four discrete velocity directions ${\boldsymbol{e}}_{1 }$ , ${\boldsymbol{e}}_{2 }$ , ${\boldsymbol{e}}_{3 }$ and ${\boldsymbol{e}}_{4 }$ is introduced and thus the lattice BGK equation for Eq (3) is written as ^[34].

where ${\tau }_{T}$ is the dimensionless relaxation time, ${T}_{i}$ is the temperature distribution function, and ${T}_{i}^{eq}$ is the equilibrium distribution function given by

The fluid temperature T can be calculated from the temperature distribution function

Similar to the recovering procedures of Eq (2), through a multi-scaling expansion, the macroscopic Eq (3) can be recovered from Eq (11) to the $O\left(\delta {t}^{2}\right)$ order if Mach number, M, is of the same order of $\delta t$ or higher, where the diffusivity $\theta$ is determined by

In the evolution of Eq (11), the thermodynamic boundary conditions with specific temperatures on the inner surfaces of the AC and PC are evaluated with the extrapolation method ^[34].

In order to coupling Eq (6) with Eq (11) in the lattice Boltzmann simulation, the Boussinesq approximation has to be performed in LBM. In our simulations, the body force is just the buoyant force, and thus the density of the body force, $\boldsymbol{F}$ , in Eq (2) can be explicitly written as

where $\boldsymbol{g}$ is the acceleration vector of gravity. Correspondingly, after absorbing the first constant part of $\boldsymbol{F}$ into the pressure term, the macroscopic Eq (2) finally becomes

Consequently, the coupling of Eqs (2) and (3), through the Boussinesq approximation, comes down to discretizing the density of the body force, $\boldsymbol{F}$ , in Eq (15) into the force term, ${F}_{i}$ , in Eq (6). By using the forcing technology proposed in ^[38], this discretization can be performed as

where i = 2 and i = 4 refer to the direction of gravity, and when i = 2 and 4, ${\alpha }_{i} = 1$ , otherwise ${\alpha }_{i} = 0$ .

With respect to the seepage flow through the TM, a relative simpler averaging method is used for solving this flow. The TM is assumed as a porous media with a finite thickness ${ h}_{TM}$ , which locates between the AC and SC. We know that the pressure in AC is just the IOP, ${p}_{IOP}$ , and that the observed pressure in the SC, ${p}_{0}$ , equals to that in blood. If the pressure in TM, ${p}_{TM}$ is supposed to vary linearly along the thickness direction, its gradient will be a constant and equal to the pressure difference cross the TM divided by the thickness ${h}_{TM}$ , e.g., $\nabla {p}_{TM} = \frac{{p}_{IOP}-{p}_{0}}{{h}_{TM}}$ . In this mentality, the coupling procedure of the incompressible Navier-Stokes flow and the Darcy's flow can be conducted as follows: firstly, we can evaluate the IOP from evolving Eq (6); secondly, the seepage velocity can be derived from the pressure gradient in TM, $\nabla {p}_{TM}$ ; and finally the calculated seepage velocity is applied to the velocity boundary condition on the AC. At this end, the coupling of the velocity field and the seepage flow through the TM has been carried out.

3.   Lattice Boltzmann dynamics of AH

In order to well understand in vivo dynamics of AH, the flow patterns and influences of various factors, such as the AH secretion rate, TM permeability and so on, should be investigated. For this matter, very few clinical and laboratorial results are available and numerical simulations may be significant methods to address these problems. Toward this end, several cases of the AH flow in both healthy (with the TM permeability given by K = 7 × 10^-15 m²) and pathological eyes (with the TM permeability given by K = 2.3 × 10^-15 m²) ^[6] situations are simulated with the present model.

3.1.   Flow of AH in healthy eye

Due to insufficient laboratorial results, the direct verification of the proposed model may encounter with difficulties, and thus an indirect validation with now available experimental data would be performed. It is known that the available experimental data for healthy eyes are much more plentiful compared to that for pathologic ones. Therefore, we will begin with the simulations for the flow of AH in healthy eye, which will be served as the basis for further researches. Let's consider a flow of AH in an anterior segment of healthy eye in standing orientation, with the geometrical dimensions listed in Table 1 and the physical parameters in Table 2.

The simulations are performed on the domain with dimensions 588 × 1921 (in lattice units). The normal production rate of AH is taken to be 2.5 μL/min, the corresponding inlet velocity is 2 × 10^-6 m/s, the temperature on surface Γ₂ is equal to 37 ℃, which is the core body temperature. In healthy eye, the normal IOP is usually 1950 Pa, and the pressure observed in SC is 1200 Pa, which is the same as that in the blood. As for the temperature, the normal temperature difference across the AC is usually 3 ℃ ^[39]. In the flow of AH, it is well known that, if the production rate is higher than the outflow rate, AH will accumulate in the AC and PC, which should lead to an increase of IOP; conversely, if the production rate is lower than the outflow rate, AH in the AC and PC will decrease, and IOP would decline. Therefore, the flow of AH will reach equilibrium with constant IOP if and only if the production rate is equal to the outflow rate. In healthy eye, when the IOP deviates from the normal value, it will automatically adjust the outflow rate through the TM, thus bringing the IOP back to the normal value gradually. To demonstrate that this regulation mechanism of the healthy eye can be simulated by using the proposed model, a test case was carried out. Initially, the IOP is set to be a deviated value from its normal one, such as 1948 Pa or 1952 Pa, and then the program evolution is performed until the system equilibrium is reached. The simulated results for the IOP and the outflow rate are drawn in Figure 2, which indicates that, with the time increase, the IOP and the outflow rate gradually converge to 1950 Pa and the inflow rate respectively. This illustrates that, due to the negative feedback of IOP, the model can always reach the equilibrium flow of AH by itself, and therefore it is convergent and robust. By changing other parameters, the simulations may reach other equilibrium flow states. Thus the proposed model for healthy eye can be applied as the basis of the further researches on pathological eye.

Thermally driven natural convection plays an important role on the dynamics of AH in the human eye ^[21]. In order to understand the importance of this heat effect, we also perform the simulations without the temperature gradient across the AC as a comparison, and the simulated results for the flow of AH with and without temperature gradient are given in Figures 3 and 4. Figure 3 (left) shows that, without temperature gradient, the flow of AH is fully generated by the AH production rate, and is rather weak with the maximum velocity of 3.03 × 10^-5 m/s and no vortex formation, which could affect the normal physiological functions of AH. With the temperature difference of 3 ℃ across the AC (Figure 3 (right)), it is clearly that the flow of AH becomes obviously strong with the maximum velocity of 7.05 × 10^-4 m/s, in which the buoyancy is the dominant driving mechanism. Furthermore, a strong recirculation vortex with its center near the center of AC is formed. To verify the results are grid-independent, the same simulations are also performed on the domain with dimensions 392 × 1281 and 784 × 2561 (in lattice units), and the maximum velocities obtained are 7.02 × 10^-4 m/s and 7.04 × 10^-4 m/s, respectively. The maximum velocities are in good agreement with the one calculated by Heys and Barocas (7 × 10^-4 m/s) ^[39].

For further investigating the recirculation of AH in the AC, which is beneficial to understand the formation of Krukenberg's Spindle, we simulate the flow of AH in both standing and up-facing eye with different temperature gradients ΔT, and the results are listed in Table 3 and drawn in Figures 6 and 7. For vortex flow, the vortex intensity is defined as

where ${\omega }_{z}$ is the angular velocity of AH, $\mathrm{\Gamma }$ is the area of AC, ${u}_{x}$ and ${u}_{y}$ are the two velocity components of AH respectively.

From Table 3 and Figure 6 (in the standing eye), we find that, with the increase of the temperature gradient, the maximum velocity and vortex intensity increase, and its center position rises slightly along the vertical direction. However, for the cases in the up-facing eye with temperature difference of 3℃ across the AC, Figure 4 shows that, driven by the temperature gradient, different from that in the standing eye, two symmetric vortices appear in the AC with equal maximum velocity of 1.05×10^-4 m/s, about 7 times smaller than that in the standing eye. Figure 7 tells us that, just like that in the standing eye, with the increase of the temperature gradient, the maximum velocity and vortex intensity also increase, and the two centers of the vortices rise along the vertical direction as well, but accompanying with their moving toward the center of AC simultaneously. Figure 5 gives the distributions of temperature in both standing and up-facing eye with temperature difference of 3 ℃ across the AC. Figure 5 (left) shows that, in the standing eye, the natural convection currents cause the posterior cornea temperature to vary, and the symmetry about the central horizontal axis of the temperature distribution is broken. Figure 5 (right) indicates that, in the up-facing eye, although the natural convection currents also cause the posterior cornea temperature to vary, the temperature distribution is symmetrical to the mid-perpendicular of AC due to the symmetry of the two vortices.

3.2.   Influence of aqueous humor secretion rate on the flow of AH

It is well-know that the IOP can be controlled by decreasing the production rate of AH. Thus the influence of AH secretion rate on IOP or the flow of AH is significant. On the basis of the simulations for healthy eye, simulations for the flow of AH with changing secretion rates are performed. With an altered secretion rate of AH, the simulation will converge to an equilibrium state with a new IOP or outflow rate. In order to understand the details of this influence on the flow of AH in both healthy (with the permeability of Κ = 7 × 10^-15 m²) and pathologic (with the permeability of Κ = 2.3 × 10^-15 m²) eyes ^[6], we investigate the two corresponding cases with changing secretion rates. The simulated results are shown in Table 4 and Figure 8. For comparisons, we also list the results obtained by ^[6] in Table 4, in which the relative error is defined as $\frac{{V}_{s}-{V}_{ref}}{{V}_{ref}}\times 100\%$ , where ${V}_{ref}$ and ${V}_{s}$ are the IOP values by ^[6] and the present model respectively.

From Table 4 and Figure 8, we see that our simulated results are very close to those obtained in ^[6], and with decrease of the inlet velocity, the IOP linearly reduces in both healthy and pathologic situations. Furthermore, the velocity of IOP reduction is quit faster in pathologic eye than in healthy eye. This numerically illustrates that the IOP can be effectively lowered by decreasing the AH secretion rate with drugs, especially in pathologic situations.

3.3.   Influence of resistance of aqueous humor drainage on the flow of AH

We also know that the IOP can be controlled by reducing the level of resistance of AH drainage. On the basis of the simulations for healthy eye, simulations for the flow of AH, with changing permeability of TM and keeping the other parameters unchanged, are carried out. When the permeability of TM has changed, the simulation will reach a new equilibrium state with the corresponding IOP or outflow rate. For comparing with the results obtained in ^[6], the simulations with the same permeability or porosity values as those used in ^[6] are performed. The simulated results with the present model, together with those from ^[6], are listed in Table 5 and drawn in Figure 9 simultaneously.

From Table 5 and Figure 9, we see that our simulated results are well agreement with those obtained in ^[6], and with the increase of the porosity (or permeability) of TM, the IOP reduces rapidly. More specifically, the velocity of IOP reduction is non-uniform, in which the IOP suffers a slight increase when the permeability K is in the interval of [7.59 × 10^-14, 5.33 × 10^-15] (m²) and but exhibits a steep gradient when K is in the interval of [3.36 × 10^-15, 5.27 × 10^-16] (m²). This numerically illustrates that the IOP can be efficiently lowered by increasing the permeability (or porosity) of TM, especially when the permeability is below the value of 3.36 × 10^-15 m².

3.4.   Influence of the viscosity of aqueous humor on the flow of AH

Viscosity of AH is a basic physical quantity associated to transport of momentum among the layers of AH, which plays an important role in the flow of AH. However, the influences of AH viscosity on the flow of AH have not been fully addressed so far. Therefore, the investigations on the influence of viscosity on the flow of AH are still needed. On the basis of the simulations for healthy eye, simulations for the flow of AH with different viscosities are carried out. When the viscosity has changed, the simulation will eventually reach an equilibrium state with a new IOP or outflow rate. We investigate two cases of both healthy and pathologic eyes. The simulated results are shown in Table 6 and Figure 10.

From Table 6 and Figure 10, we see that, with the increase of the viscosity, the IOPs in both healthy and pathologic eyes will all rise at different velocities, and the IOPs rise faster in the pathologic eye than in the healthy one. But, on the other hand, the IOP increase is insensible to the viscosity of AH, compared to the permeability of TM and AH secretion rate.

3.5.   Influence of anterior chamber angle on the flow of AH

Anterior chamber angle is defined as the attachment angle of the iris to the cornea, which is closely related with close-angle glaucoma. In order to research the IOP dependence on anterior chamber angle, we perform the simulations with different anterior chamber angles in both healthy and pathologic eyes, in which the range of the angle is much larger than that in ^[15]. The simulated results are listed in Table 7 and drawn in Figure 11. Both Table 7 and Figure 11 indicate that, with the decrease of the anterior chamber angle, the IOPs increase in both healthy and pathologic eyes with increasing velocities, of which the IOP increase velocity is slightly greater in the pathology eye than in the healthy one. Analyzing in detail, we find again that the IOP increases slowly with the decrease of the anterior chamber angle when this angle is above some critical value (about 15 ℃), and it dramatically increases once the angle is lower than the critical value, which implies a higher risk of glaucoma.

3.6.   Wall shear stress

The flow of AH in the anterior chamber will form a wall shear stress (WSS) on the corneal endothelium or the iris surface. Researches show that WSS can change the endothelial cells' functions and structures ^[40,41]. In this subsection, the factors that may affect these wall shear stresses are investigated in a healthy eye. The shear stress contour for the standing and up-facing orientations of the eye are plotted in Figure 12 with a temperature difference across the AC ΔT = 3 ℃. It can be seen that there is an obvious difference for the distribution of WSS on the corneal endothelium and the iris surface between the standing and the up-facing orientation of the eye. For the standing orientation, the maximum cornea WSS is located at the top of cornea and decreases towards the anterior chamber angle, while the maximum iris WSS is located in the inner periphery of the iris close to the pupil and lessens towards the anterior chamber angle. For the up-facing orientation, the maximum cornea WSS is located in the midperiphery of the cornea and reduces towards the anterior chamber angle and the top of cornea respectively, and the maximum iris WSS appears in the middle of the iris surface and decreases towards the pupil and the anterior chamber angle separately.

Table 8 shows the dependence of average cornea and iris WSSes with the temperature difference ΔT. It can be seen that there is a linear relationship between the WSSes and the temperature difference. The greater the temperature difference, the greater WSS is. The WSS is sensitive to the temperature difference, when the temperature difference increases from 0.2 ℃ to 7 ℃, there is an increasing more than an order of magnitude on the cornea or iris WSS in both standing and up-facing orientations of the eye. For instance, the cornea WSS increased by about 45 times in the up-facing orientation. Table 9 shows the dependence of average cornea and iris WSSes with the inflow velocity (AH secretion rate), permeability of TM and AH viscosity. For the inflow velocity, the average WSSes in standing and up-facing orientation increase by 0.05% and 1.75% for the cornea, 2.19% and 4.16% for the iris respectively when it increases by 4 times. The increase rate of the average WSS in the up-facing orientation is slightly higher than that in the standing orientation. For the permeability of TM, according to Eq (4), it is directly proportional to the outflow velocity of AH. However, we found that the average cornea and iris WSSes remained almost unchanged when the permeability of TM increases by 144 times. The cause may be that the gradient of AH velocity near the wall of the cornea and iris do not increase although the outflow velocity of AH increases. In order to investigate the effect of the AH viscosity on WSS, we increased the AH viscosity by 1.23 times, the average cornea and iris WSSes in the standing orientation increase by 1.23% and 3.69%, while those in the up-facing orientation decrease by -1.14% and -1.40%, respectively. According to Newton's law of viscosity, the viscosity and velocity gradient of AH are in direct proportion to the shear stress. However, the velocity gradient of AH near the wall decreases when the viscosity of AH increases, so the average cornea and iris WSSes in the standing orientation have a slight increase. On the contrary, the WSSes in the up-facing orientation have a little reduce due to the influence of the velocity gradient decreasing on the WSS is greater than that of the AH viscosity. The same simulations mentioned above were also performed in a pathologic eye. However, we found that the results are almost the same as those in the healthy eye. The cause is that the healthy and the pathologic eyes are defined by the permeability of TM K = 7 × 10^-15 m² and K = 2.3 × 10^-15 m² respectively in this paper, and the above results have proved that the permeability of TM has little effect on the WSSes. Through analyzing above simulation results, we can draw conclusions that the average cornea and iris WSSes in the standing orientation is significantly higher than those in the up-facing orientation, and among the factors affecting the WSS, the temperature difference has the greatest influence on the WSSes, while the influence of inflow velocity, permeability of TM and AH viscosity are not obvious.

The magnitude of the WSSes on the iris surface and corneal endothelium is an important parameter for the prevision of endothelial cell detachment in pigment dispersion syndrome, pigmentary glaucoma and bullous keratopathy ^[42,43]. In our simulations, the average the maximum cornea and iris WSSes are 1.61 × 10^-3 Pa and 1.85 × 10^-3 Pa respectively. According to Gerlach et al. ^[44], the magnitudes required to detach endothelial cell on the iris surface are in the rage between 0.51 and 1.53 Pa, while the ones required to detach endothelial cell on the corneal endothelium are comprised between 0.01 and 1.0 Pa as reported by Kaji et al. ^[45]. Due to the WSSes in our cases are much lower than the detaching threshold value for such cells, it is unlikely that the WSS could detach the pigments from the anterior iris surface or the endothelial cells from the corneal endothelium, but the plots of shear stress give an insight about the location and probability of pigments detachment on the iris surface or corneal endothelial cells loss on the corneal endothelium in case of some specific disease, such as angle-closure glaucoma after laser iridotomy.

4.   Conclusions

A coupled lattice Boltzmann model for simulating the dynamics of aqueous humor in the human eye is presented, which involves incompressible Navier-Stokes flow, heat convection and diffusion, and Darcy seepage flow. The verification simulations for the healthy eye show that the present model is available, convergent and robust due to its pressure negative feedback mechanism. Further researches on the dynamics of AH in the healthy eye indicate that heat convection and diffusion is crucial to the normal physiological functions of AH. Furthermore, the flow patterns of AH, driven by the temperature differences across the AC, are quite different between the standing and up-facing eyes. The temperature differences and the orientations of the eye have significantly influences on the flow patterns of AH, including the maximum velocity, the numbers of the vortex, the intensities and center positions of the vortex and so on, which may be beneficial to understand the formation of Krukenberg's Spindle.

Basing on the verification simulations, we systematically research the influences of the AH secretion rate, permeability of TM, viscosity of AH and anterior chamber angle on IOP in both healthy and pathologic eyes, the WSSes on the corneal endothelium and the iris surface, and the following conclusions can be drawn from the detailed analysis of the results.

(1) With the decrease of the inlet velocity or secretion rate of AH, the IOP linearly reduces in both healthy and pathologic eyes. Furthermore, the IOP reduction exhibits a steeper gradient in pathologic eye than in healthy one. Therefore, through reducing the AH secretion rate with drugs, the IOP may be expected to be controlled efficiently.

(2) With the decrease of the TM permeability, the IOP will gradually increase with increasing velocity. In more detail, there would be a critical value of the TM permeability (about 5.33 × 10^-15 m²), above which the IOP increases slowly and below which the IOP will increase rapidly. This finding is consistent with the clinical practice, in which the diagnosis of glaucoma does not simply rely on the presence of high IOP.

(3) With the increase of the AH viscosity, the IOPs in both healthy and pathologic eyes will all rise at different velocities, and the IOP rises faster in the pathologic eye than in the healthy one. However, the IOP increase is insensible to the viscosity of AH, compared to the permeability of TM and AH secretion rate, which may be one of the causes that the viscosity would not have been considered as a factor for controlling IOP.

(4) With the decrease of the anterior chamber angle, the IOPs increase in both healthy and pathologic eyes with increasing velocities. For the IOP increase rate, there would be some critical value of the anterior chamber angle (about 15 ℃), above which the IOP increases slowly and below which the IOP increases dramatically. In addition, the IOP increase is faster in the pathologic eye than in the healthy one, which indicates that the narrowing of the anterior chamber angle produces a more significant effect in a pathologic situation.

(5) The WSSes in the standing orientation is notable higher than those in the up-facing orientation, and the temperature difference have significant influence on the cornea and iris WSSes, while the influence of inflow velocity, permeability of TM and AH viscosity are not obvious. The effect on WSSes with all these factors in pathologic eye is almost same as those in healthy eye. The WSSes in our cases are much lower than the detaching threshold value for endothelial cell, so it is unlikely that these WSSes could detach pigments from the anterior iris surface or endothelial cells from the corneal endothelium.

Due to the advantages of the lattice Boltzmann method in modeling complex fluid systems, the present model can be expected to become a good alternative tool for studying the AH dynamics and apply to the dynamics of AH more deeply and obtain some more interesting results.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant Nos. 12062005, 81860635, 11862003), the Project of Guangxi Natural Science Foundation (Grant Nos.2017GXNSFDA198038, 2018GXNSFAA281302), the Graduate Innovation Program of Guangxi Normal University (Grant No. XYCSZ2020070), and the Guangxi Collaborative Innovation Center of Multi-source Information Integration and Intelligent Processing.

Conflict of interest

All authors declare no conflicts of interest in this paper.