- Research
- Open Access
- Published:

# Development of a 3D finite element model of lens microcirculation

*BioMedical Engineering OnLine*
**volume 11**, Article number: 69 (2012)

## Abstract

### Background

It has been proposed that in the absence of a blood supply, the ocular lens operates an internal microcirculation system. This system delivers nutrients, removes waste products and maintains ionic homeostasis in the lens. The microcirculation is generated by spatial differences in membrane transport properties; and previously has been modelled by an equivalent electrical circuit and solved analytically. While effective, this approach did not fully account for all the anatomical and functional complexities of the lens. To encapsulate these complexities we have created a 3D finite element computer model of the lens.

### Methods

Initially, we created an anatomically-correct representative mesh of the lens. We then implemented the Stokes and advective Nernst-Plank equations, in order to model the water and ion fluxes respectively. Next we complemented the model with experimentally-measured surface ionic concentrations as boundary conditions and solved it.

### Results

Our model calculated the standing ionic concentrations and electrical potential gradients in the lens. Furthermore, it generated vector maps of intra- and extracellular space ion and water fluxes that are proposed to circulate throughout the lens. These fields have only been measured on the surface of the lens and our calculations are the first 3D representation of their direction and magnitude in the lens.

### Conclusion

Values for steady state standing fields for concentration and electrical potential plus ionic and fluid fluxes calculated by our model exhibited broad agreement with observed experimental values. Our model of lens function represents a platform to integrate new experimental data as they emerge and assist us to understand how the integrated structure and function of the lens contributes to the maintenance of its transparency.

## Background

Our sense of sight is critically dependant on the ability of the optical pathway, formed by the cornea and lens, to focus light onto the retina. As an optical element in this pathway, the lens needs to maintain its transparency and create sufficient optical power. To generate its required optical properties, the lens has evolved a unique structure to minimise light scattering and enhance transparency [1] [see Figure 1. However, the lens is not a purely passive optical element and its structure and therefore optical properties need to be actively maintained. It has been proposed by Mathias et al. [2] that in the absence of a blood supply, the lens operates an internal microcirculation system. This system delivers nutrients, removes wastes and maintains the ionic homeostasis of the lens [3, 4]. Hence, this system is critical for maintaining the optical properties of the lens.

The microcirculation model was initially based on a combination of electrical impedance measurements and theoretical modelling
[8–11]. The observation that the measured ionic currents were directed inward at the poles and outward around the equator [Figure
1A, led to the suggestion that these currents represent the external portion of a circulating ionic current that drives the internal microcirculatory system within the lens. Briefly, the working model is that the current, carried primarily by Na^{+}, enters at all locations around the lens via the extracellular space between fibre cells. Na^{+} eventually crosses the fibre cell membranes, and then flows from cell-to-cell towards the surface via an intracellular pathway mediated by gap junction channels [Figure
1B. The gap junction coupling conductance in the outer shell of differentiating fibres is concentrated at the equator
[12, 13]. Hence, the intracellular current is directed to the equatorial epithelial cells. Here, the highest densities of Na^{+}/K^{+} pumps are located to actively transport Na^{+} out of the lens
[14]. Thus, the intracellular current effluxes are highly concentrated at the equator causing the net current to be outward. At the poles, there is very little intracellular current so the net current is predominantly inward, along the extracellular spaces [Figure
1B. The driving force for these fluxes is hypothesized to be the difference in the electromotive potential of surface cells that contain Na^{+}/K^{+} pumps and K^{+}-channels, and inner fiber cells that lack functional Na^{+}/K^{+} pumps and K^{+}-channels and whose permeability is dominated by non-selective cation and Cl^{-} conductances
[15]. This electrical connection together with the different membrane properties of the surface and inner cells causes the standing current to flow. In this model, the circulating current creates a net flux of solute that in turn generates fluid flow. The extracellular flow of water convects nutrients towards the deeper lying fiber cells, while the intracellular flow removes wastes and creates a well-stirred intracellular compartment.

Thus in this model, it is the circulating Na^{+} current that generates a circulation of fluid inside the lens. It is important to note that while the existence of the circulating ionic currents are firmly based on existing experimental data, circulating fluid flows in the lens have proven more difficult to measure. At present, these fluxes are only predicted to occur from indirect measurements and models of the measured electrical properties of the lens. The current model of the lens fluid dynamics is based on an equivalent circuit analysis of the microcirculation system
[3, 16]. This analytical model inherently relied on approximate solutions and simplification of the underlying physics. To improve our understanding of the circulation system, here we have adopted a finite element modelling (FEM) approach. Using FEM, we have developed a 3D computer model of the microcirculation system that encapsulates the complex interplay of its important features and can also be solved numerically. This computer-based modelling approach allowed structural features such as fiber cell orientation, extracellular space dimensions and gap junction distribution to be included in the model. It also included functional information on the spatial differences in membrane permeability between surface and inner cells, thought to drive the circulating currents.

In this paper we describe our expansion of the equations that govern ion and fluid dynamics in the lens [2, 3, 10, 15–18] to 3D, and their subsequent implementation into a new FEM that encapsulates known structural and functional parameters of the mouse lens. This model is based on our continuous imaging and modelling iterative investigation into the fluid dynamics of the lens [19–25]. To test the ability of this computer model to reproduce the functional properties of the lens, we have used a series of experimentally derived boundary conditions [26–28] to allow it to be solved. We show that our model is capable of predicting the experimentally measured steady state lens properties and to generate circulating ion and water fluxes. Hence we believe that our computer model is a useful tool to study how lens structure and function influence its optical properties.

## Methods

In this section, we first present the derivation of the fundamental mathematical equations, originally formulated by Mathias et al. to describe the microcirculation system [2, 3, 10, 15–18]. We then develop a computer mesh to represent the structure of the mouse lens to enable these equations to be implemented in 3D. Next we implement the model using the C++ programming language and solve the model using experimentally derived boundary conditions [26–28]. The resultant 3D model calculates ion concentration gradients, membrane potential gradients, and intra- and extracellular ion and water fluxes in different regions of the lens.

### Derivations of general equations

The microcirculation model has been based on a distinct yet mutually dependent movement of water and ions [3, 4, 15, 16, 29]. Simplified derivations of the Navier–Stokes equations and the advective Nernst-Plank equations [30, 31] where used to represent water and ion fluxes, respectively. All symbols used in these and the resultant equations are listed in [Table 1.

#### Fluid fluxes

The Navier–Stokes equations are derived from the conservation of mass, momentum, and energy. The nonlinearity of these equations made them challenging to be solved directly. However, this set of equations can be simplified with the right assumptions for a given system. We assumed the lens’s water to be an incompressible Newtonian fluid with a spatially constant viscosity at steady state. Also, we assumed that the lens fluid flow can be described as a “creeping” low-Reynolds number flow with ignorable turbulence [32–34]. Considering these assumptions, we used the simplified form of the general Navier–Stokes equations, called the Stokes equations [Eq. 1, Eq. 2.

The parameters and their units are listed in [Table 1. We used these equations to model the movement of water in the lens. All the fluid flow related constants utilized in this model are listed in [Table 2. Most of the parameters are similar to the parameters used in the analytic model [2, 3, 10, 15–18]. The most significant change between the analytic and numerical models is the intercellular, surface and fibre cell hydraulic permeabilities, where more recent values, which were measured by [35], are used.

#### Ion fluxes

The solute fluxes in a fluid are generally governed by diffusion, electro-diffusion (if the solute is charged) and advection [30, 31, 37]. Hence, we modelled the ionic fluxes in the lens using the Nernst-Plank equation with an added advection term [Eq. 3 [18, 30, 31, 37].

The parameters and their units are listed in [Table
1. All the ionic flow related constants utilized in this model are listed in [Table
3. The intracellular diffusion coefficients for Na^{+}, K^{+} and Cl^{-} (D_{i,c}) in the radial direction are estimated to be 1% of the diffusion coefficients in the cytoplasm
[38]. The diffusion coefficient for small ions in the cytoplasm is assumed to be the same as in free solution. This assumption is valid if there is no significant interaction between the ions and the cell membrane or large molecules. Similarly, the extracellular diffusion coefficients for Na^{+}, K^{+} and Cl^{-} (D_{e,c}) are assumed to be the same as in free solution. These values are scaled by the τ_{c} [Table
3 to account for the tortuosity of the extracellular cleft. The conductivities for the Na^{+}, K^{+} and Cl^{-} channels are the same as those used in the analytic model and are assumed to be spatially uniform.

We solved these equations [Eq. 1 to Eq. 3] for the movement of ions and water in the extracellular and intracellular spaces of the lens. These two domains were linked by cell membranes which could be crossed by ions and water from one space to another [Figure 1B]. Since the water and ions enter the lens from the extracellular space, we started with these equations.

### Extracellular flux equations

The fluid flow in the extracellular clefts can be partially described by the Stokes flow equations [Eq. 1 & Eq. 2. Another part of the extracellular fluid fluxes has been shown to be due to electro-osmosis [39]. Electro-osmosis is due to the osmotic gradient created by uneven charge distribution in a fluid affected by an electric field. It has been shown that this osmosis is essential to the modelling of extracellular fluxes in the ocular lens [39].

*Extracellular fluid flows* in the lens were modelled by the Stokes equations [Eq. 1 & Eq. 2 with an added electro-osmosis term. We treated electro-osmosis as a body force term and formulated it as below for our lens model [Eq. 4
[18].

The parameters and their units are listed in [Table
1]. We then combined this velocity (*u*
_{
eo
}) and the Stokes equation [Eq. 2] to drive the equation for the extracellular water velocity [Eq. 5].

The parameters and their units are listed in [Table 1. We didn’t consider electro-osmosis in the intracellular space portion of our model [18]. This was since the electric field gradient across cell cytoplasm was considered negligible. We then coupled the extracellular water fluxes with the ionic fluxes in this domain.

#### Extracellular ion fluxes

We computationally solved the advective Nernst-Plank equation [Eq. 3] for the lenticular extracellular ion fluxes. We modelled the coupled water and ion extracellular fluxes throughout the 3D model. These fluxes could then become trans-membranous flows at any point of the model, crossing the modelled cell membranes.

### Trans-membrane flux equations

The extracellular and intracellular spaces of this model are connected by cell membranes with specific water and ionic permeabilites. Hence the water could flow cross fibre cell membrane between the extracellular and the intracellular spaces. To represent these transmembrane fluxes in our model, we treated the fibre cell membrane as a semi-permeable membrane [2, 7, 10, 16] through which fluid passes due to a combination of hydrostatic and osmotic pressure gradients [19, 40, 41]. The velocity of the water flowing through the membrane was then governed by the following equation [Eq. 6 [18].

The parameters and their units are listed in [Table
1]. We calculated the hydraulic and osmotic pressure gradients, *Δp* and *ΔOs*, from the pressure differences on two sides of the membrane. We estimated the general osmolarity, *Os*, of the intracellular or extracellular spaces from the local sum of all the modelled ions (i.e. Na^{+}, K^{+} and Cl^{-}).

As mentioned previously, we modelled the ion and water fluxes to go hand-in-hand throughout our computational platform. Hence, the modelled ions (i.e. Na^{+}, K^{+} and Cl^{-}) accompanied the trans-membrane water fluxes into the cells, utilizing the membrane ionic channels. The lenticular membrane conductivity for each modelled ion had been calculated based on experimental data
[8, 11, 42, 43]. We used similar values for membrane conductivities for various modelled trans-membrane ion fluxes. We implemented these trans-membranous ionic flows using the following equations [Eq. 8 to Eq. 10
[18].

The parameters and their units are listed in [Table
1]. In the above equations, *E* is the Nernst potential. At the Nernst potential, there is no net flow of ions across the cell membrane through the channels. Hence, in our model we linked the trans-membrane water and ion fluxes by ionic concentrations (i.e. osmosis) and membrane potential (i.e. Nernst potential). After crossing the membrane, we treated the water and ion fluxes as intracellular fluxes.

### Intracellular flux equations

We modelled the fluxes in the intracellular space to pass through a mixture of cells cytoplasm and membranes [44]. This was due to the large size of each element in our model compared to the lens cell volumes. It was impractical to discretely model the intracellular flow through each cell cytoplasm and membrane in a given element. Instead, we homogenized the flow equations to obtain one formula, describing the net cytoplasmic and membranous fluxes through every element.

To homogenize the flow equations, we assumed the intracellular fluid fluxes are primarily restricted by the cell membranes. We made this assumption since it was shown that diffusivity between cells is about on hundred times lower than in the cytoplasm [44]. We also assumed negligible hydraulic and osmotic pressure gradients across the cytoplasm of a single cell. We considered these fields to vary only among neighbouring elements, across the modelled membranes. Hence we assumed the pressure fields as step-function changes across the cell membranes. We therefore implemented the intracellular fluid flow to be driven by the pressure and concentration step-changes across the cell membranes. We calculated the intracellular fluid velocities using the following equation [Eq. 11.

The parameters and their units are listed in [Table
1. Since the implemented intracellular water fluxes were accompanied by the ionic fluxes in our model, we assumed the intracellular ionic fluxes to be governed by the advective-Nernst-Plank equation [Eq. 3. As mentioned above, every element of the intracellular space in our model was treated as a homogenised mixture of cytoplasm and membranes. In reality a network of gap junctions forms cell-to-cell channels that connect the cytoplasms of adjacent cells
[40, 45]. We derived the ionic flux equation for a single gap junction channel [Eq. 3 and then integrated this equation over a network of pores in a given cell membrane
[18]. We used [Eq. 3 for ionic fluxes, assuming free movement of ions through each pore. In this case we could assume the diffusion coefficient (*D*) in the pore to be equivalent to the cell’s interior diffusion coefficient. We also assumed the ion concentrations (*C*) and electric potentials (*ϕ*) to vary linearly between the opposite openings of each pore. Therefore, we could approximate these variables over the length of the pore (*x*), by the following equations [Eq. 12 & Eq. 13.

Here, the subscripts *1* and *2* referred to the inlet and outlet of the modelled pore. We approximated the ion concentration in the pore (*C*
_{
α
}) by the mean concentration [Eq. 14].

By substituting the above equations [Eq. 12 - Eq. 14] in the advective-Nernst-Plank formula [Eq. 3], we derived the following form [Eq. 15].

The subscript *α* indicated the solute species and the rest of the variables are listed in [Table
1]. The diffusion coefficient (*D*) in the above equation is mainly set by the directional conductivity of the gap junctions [Eq. 16].

The parameters and their units are listed in [Table
1. In this implementation, the diffusivity of a gap junction pore (*D*
_{
α
}) is related to its conductivity (*G*
_{
α
}). Based on experimental observations, the measured conductivities of the mouse lens have been assorted into two main regions. The inner 85% of the lens is anucleated
[1] and has an approximately spatially uniform distribution of gap junctions. For this inner part of the lens a homogeneous value of conductivity (*G* = 0.17mS/cm) has been used
[7]. The outer 15% of the lens on the other hand is nucleated
[1]. The radial conductivity, G_{r} has been assumed to vary depending on the angular location (*ϕ*) and the conductivity of the gap junctions [Eq. 17
[7].

Here *ϕ* is the angular degree from the equator (i.e. *ϕ* = 0° at equator, *ϕ* = 90° at anterior pole and *ϕ* = −90° at posterior pole). It has been determined that *G*
_{
max
} = 0.6mS/cm in the radial orientation, and in the other directions *G* = 0.23mS/cm
[7]. We modelled the intracellular ionic fluxes to be directed by the above equations and the regional distribution of gap junctions [Eq. 15 - Eq. 17. According to the microcirculation theory, these ionic fluxes then moved towards the periphery of the lens accompanied by water flows until they reach the boundary of the model.

### Surface flux equations

Once intracellular fluxes reach the boundary of our 3D computational mesh they are treated as surface fluxes. We modelled the periphery of the lens to be similar to previously implemented cell membranes. The only difference was that at this point we had to consider the surrounding boundary conditions that are determined by the ionic composition of the bathing media in the system of equations. At the surface of the lens, fluid is transported between the intracellular space and the outside. We modelled these surface fluid flows as simple trans-membrane fluxes and calculated their velocities using the equation [Eq. 18].

The parameters and their units are listed in [Table
1] and the subscript “*o”* points to the outside (i.e. boundary) conditions. We considered the surface fluxes as the output of the system, while the extracellular fluxes were its input. For an incompressible stable system, like the current lens model, the input and output levels should equate at all times and this conservation of mass was controlled for as part of the Stokes equations [Eq. 1].

At the surface of the lens, we modelled the intracellular space to be coupled to the outside media. We created this coupling by a network of K^{+} channels and Na^{+}/K^{+}-pumps in the surface membranes
[39, 44]. To model ionic fluxes through the K^{+} channel, equations [Eq. 8 to Eq. 10 were used. To implement the Na^{+}/K^{+}-pumps at the surface of the lens we used a previously published approach
[39] that utilised the following set of equations [Eq. 19 to Eq. 21 to model these pump rate.

Here, *I*
_{
max
} was the maximum current rates through the pumps. The parameters and their units are listed in [Table
1. The *K*
_{
0.5Na
} and *K*
_{
0.5K
} were the concentrations of Na^{+} and K^{+} when the pump’s current was half the maximum current
[39]. We also implemented the known circumferentially-varying distribution of ionic pumps around the lens
[7]. The ionic pumps of the lens are more abundant around its equatorial area and scarce around the polar regions
[7, 39]. It has been suggested that the currents produced by Na^{+}/K^{+} ATPase pumps around the ocular lens are governed by [Eq. 22
[7].

Here *ϕ* is the angular degree from the equator (i.e. *ϕ* = 0° at equator, *ϕ* = 90° at anterior pole and *ϕ* = −90° at posterior pole). Experiments have shown that *I*
_{
p-min
} = 4 μA/cm^{2}
[39] around the poles and *I*
_{
p-max
} = 26 μA/cm^{2} near the equator
[41]. Hence, at this point the model included the regional distribution of the ionic pumps. We then calculated the surface Na^{+} and K^{+} fluxes by using the pump rate (*I*
_{
p
}) estimated above. We estimated these pump-generated molar fluxes at every point on the surface of the 3D mesh [Eq. 23 & Eq. 24.

We modelled the total surface ionic fluxes by combining the pump-generated and channel-based flows [Eq. 8, Eq. 23 & Eq. 24].

### Electro-neutrality equation

Being a closed system, the model operates in the condition that at every point of time electro-neutrality is preserved. In this model we used a weak form of electro-neutrality [Eq. 25] to enforce this constraint.

We assumed the initial condition of the system to be electrically neutral. Hence, the equation above ensures the electro-neutrality of the model at any point of time. We applied this equation over all the ion species modelled here. This was based on the assumption that the ion species not modelled (e.g. calcium) had no significant influence on electro-neutrality of the model.

### Solving the model

We implemented all of the mentioned water and ion flux equations in our model. This was in order to create an interlinked system of fluid dynamics of the ocular lens. The implemented system of equations is summarized in the following figure [Figure
1B. To solve this system of equations we used an adaptive Euler method, used to achieve a converged steady state solution
[18]. Each iteration of the adaptive Euler method involved several steps to solve the coupled transport equations, followed by a solution update step [Figure
2. The model began with a representative mesh of the lens, the initial conditions (*C*
_{
0
}
*, ϕ*
_{
0
}) and boundary conditions (*C*
_{
αo
}
*, ϕ*
_{
o
}
*, p*
_{
o
}). Ideally, one would start with initial concentration fields that are close to the expected final concentration fields.

### Finite Element Mesh Creation

Since a wealth of experimental data exists on the mouse lens we chose it to create an anatomically accurate scaffold to implement our modelling approach. The mouse lens has an equatorial radius of 0.125 cm, posterior thickness of 0.1 cm and anterior thickness of 0.085 cm [42]. We used a cylindrical polar coordinate system (r, θ, z) and the Cubic Hermite basis function to create a smooth 3D computational mesh of the mouse lens [Figure 3 and Table 4. When a biological tissue is modelled as a macroscopic continuum, elements in the representing mesh are much larger than the cells. In our model, the intra- and extracellular spaces of the lens therefore could not be represented disjointedly. Hence, we modelled that these spaces, while mathematically distinct, coexist in the same mesh element and the cell membranes are evenly spread all through. Every element in our model thus represented a cluster of many fibre cells and enclosed extracellular space. This method of “bidomain” modelling has been well-explained previously [43].

### Boundary conditions

Since our model is solved using an iterative method, it was important to choose appropriate initial boundary conditions. In previous microcirculation models, the boundary conditions have been chosen to be ionic concentration at the surface of the lens (*C*
_{
Na-io
}
*, C*
_{
K-io
}
*, C*
_{
Cl-io
}
*, C*
_{
Na-eo
}
*, C*
_{
K-eo
}
*, C*
_{
Cl-eo
})
[3, 10]. Here subscript “o” denotes the surface of the lens. We adopted the same approach to solve our 3D model. Boundary conditions were taken from the measurements of ionic concentrations from the surface of ocular lenses in different species
[26–28, 47] and are listed in [Table
5. Using these initial conditions we solved the model and assumed the initial extracellular electric-potential and hydraulic pressures to be zero. These fields were then calculated during the first iteration round and substituted in the model for the following iteration cycle with this iteration being repeated until convergence was reached.

### Convergence criterion

We defined the maximum change of concentrations between two consecutive iterations (*C*
_{
n+1
}
*- C*
_{
n
}) in all the elements as the convergence parameter. We observed that using the above boundary and initial conditions, this error was decreased with each iteration cycle. We set the model’s convergence criterion to be less than 5 mM. We believed that level of iteration error was adequate, since most of the solution fluctuations were caused by [*K*^{+}]_{
i
} and [*Na*^{+}]_{
e
} fields and the enforced iteration error threshold was less than 5% of the initial field value.

We solved the 3D model on our high performance computer (HPC) at the Auckland Bioengineering Institute (ABI). Our mainframe was comprised of an IBM server with 64 processors of 1.9 GHz calculation speed, peak performance of 306.21Gfps and 256 GB of physical memory.

### Displaying the data

Text formatted file of the model’s computed fields were linked via JAVA programming language format to CMGUI ( http://www.cmgui.org) [48, 49], an advanced 3D visualisation software package with modelling capabilities. CMGUI was used for model visualisation and manipulation and allowed for automated the scaling and false colouring of the data. The following presentations were all created using this method.

## Results

Using the boundary conditions listed in [Table 5], we solved the model and generated 3D maps of standing fields of intracellular and extracellular ion concentrations, electrical potentials and pressure, plus circulating ionic and water fluxes. In this section we first use quarter section views of the 3D model [Figure 3C] to represent regional differences in standing electrochemical and pressure fields, allowing these predicted properties to be compared to experimentally derived values. Then we use the full model view to generate 3D vector maps that visualize the predicted ion and water fluxes in the lens for the first time.

### Standing electrochemical and pressure fields

The operation of the microcirculation system will at steady state create standing fields in ion concentrations and membrane potentials in both the extracellular and intracellular compartments. The 3D representations of the ion concentration fields for Na^{+}, K^{+} and Cl^{-} in both compartments are illustrated in Figure
4. The values for these concentrations were then extracted and presented in Table
6. To assess the accuracy of the standing ionic gradients produced by the model, we wanted to compare them to existing experimental data. Unfortunately no data on the concentration of ions in the extracellular space of the lens is available. However, ion concentrations in the intracellular space have been measured using a variety of techniques. Wang et al.
[50] have used ion selective micro-electrodes to show a Na^{+} concentration gradient exists in the mouse lens, the magnitude of which is predicted by our model [Table
6. Average concentrations of Na^{+}, K^{+} and Cl^{-} in the lens have also been measured using chemical analysis in a number of species
[26–28, 47]. To facilitate comparison between these previous measurements of total lens ion concentrations and the values predicted by our model, we extracted the mean concentration from our calculated ion fields [Table
6. Here we observed a good compatibility between the modelled and measured values.

The trans-membrane ionic concentration fields determined by the model were then used to calculate the intra- and extracellular and trans-membrane potential fields [Figure
5. These calculated values agreed with lens potentials measured experimentally
[9, 51], [Table
6. The model predicted a decrease in extracellular space electrical potential from the periphery to the core. This is thought to be indicative of an increase of extracellular space resistance with depth into the lens
[8]. In contrast, a smaller change in the intracellular electric potential field was predicted by the model. This small change has been observed previously
[2, 10] and has been attributed to a high level of gap junction coupling between inner-cells and surface cells containing K^{+} channel and pumps that dominate the lens potential
[15]. This meant that the potential between the extracellular and intracellular spaces depolarized from -63 mV at the surface to -33 mV in the centre of the model [Figure
5C. These changes in potential fields with distance into the lens were consistent with the movement of a positively-charged ion like Na^{+}, into the lens via the extracellular space, down its electrochemical gradient across the cell membrane and outwards in the intracellular space.

The Stokes equations [Eq. 1 & Eq. 2 used to model the fluid dynamics of the lens has an inherent pressure term, which allows us to calculate the standing fields of hydrostatic pressure within the intra and extracellular compartments [Figure 6. The model calculated an extracellular pressure gradient that was highest at the periphery and lowest in the core [Figure 6A and an intracellular hydrostatic pressure gradient that was lowest in the periphery and highest in the core [Figure 6B. Recently the existence of an intracellular pressure gradient has been confirmed experimentally [29]. A comparison of the calculated and measured intracellular pressure fields indicated that the orientations of the gradients are similar but the magnitudes differ by a factor of ~4.5 [Table 6. This disagreement is discussed further in the Discussion.

The results obtained from our 3D model of the calculated electrochemical and pressure fields can be simplified to facilitate comparison between parameters, by extracting values from a defined equatorial axis and plotting them against normalized distance from the centre (r/a) [Figure 7]. In the lens, these fields give rise to the circulating fluxes modelled below.

### Circulating fluxes

We calculated extracellular and intracellular ionic fluxes throughout the lens using the advective Nernst-Plank equation. Extracellular space fluxes of Na^{+} were directed into the lens and had the highest magnitude at the poles [Figure
8A
[4]. In contrast, intracellular space fluxes of Na^{+} were directed out of the lens and had their highest magnitude at the equator [Figure
8B. Extracellular [Figure
8C and intracellular [Figure
8D fluxes of Cl^{-} appeared to mirror the Na^{+} fluxes as would be expected to preserve electro-neutrality. Relative to extracellular fluxes of Na^{+} and Cl^{-}, the calculated extracellular K^{+} fluxes were of reduced magnitude [Figure
8E. However the intracellular space K^{+} fluxes were of larger magnitude and were concentrated at the equator [Figure
8F, consistent with the abundance of pumps and K^{+} channels known to be localized in this region of the lens
[3, 10, 16].

Although ionic fluxes within the lens have yet to be determined experimentally, measurements of net current densities at the lens surface are available [5, 41]. To facilitate a comparison between these measurements and our model’s predictions, we calculated net current densities using [Eq. 26 & Eq. 27.

The parameters and their units are listed in [Table 1. The resultant 3D map of the current densities throughout the lens are shown in Figure 9A, and the values extracted for the surface densities located at both poles and the equator are listed in Table 7. Our model clearly shows a net current influx at the poles and a net efflux at the equator [Figure 9B, a pattern in agreement both in direction and magnitude with previous electrophysiological measurements conducted on a variety of lenses [Table 7, [5, 41]. What our model also showed for the first time was the pattern in magnitude of net current fluxes within the lens.

The calculated ionic current densities however, are only one leg of the microcirculation system. To study the water fluxes generated by ion flow, we used the Stokes equations to calculate water velocities in the intracellular and extracellular spaces of our model [Figure 10. We observed that the extracellular fluid velocities were inward and maximal at the poles [Figure 10A. Conversely, the intracellular fluid velocities were outward and maximal at the equator [Figure 10B. Thus our model generates the theoretical water fluxes first proposed by Mathias et al. [16], but unfortunately no experimental data currently exists to verify the accuracy of the calculated magnitudes [Table 7.

## Discussion

It has been proposed that the circulating currents observed experimentally in lenses from a variety of species [5, 14, 41, 52] drive an internal microcirculation system. In the absence of a blood supply the microcirculation delivers nutrients and removes metabolic waste products from inner fibre cells, while maintaining ionic homeostasis [7]. The circulating fluxes are thought to be the net result of spatial differences in the location of ion channels and transporters that determine local membrane permeability, and the density of gap junction channels that direct intercellular fluxes within the lens [3]. We first created a 3D mesh of the lens structure. We then implemented a series of equations that describe the local transport properties of cells in different regions in the lens. Finally we solved these equations using a FEM approach at each location of the model. In summary, we have produced a computer model of the lens that not only accurately predicts standing ionic concentrations [Figure 4 and electrical [Figure 5 gradients and the existence of a pressure gradient [Figure 6, but also produces the first 3D vector maps of predicted current [Figure 8 and fluid [Figure 10 fluxes inside the lens. The observed agreement between experimentally measured values and those calculated by our model [Table 6 & Table 7, suggest that our computer model is mimicking lens physiology and generates a microcirculation. However, the underestimation of the magnitude of the intracellular pressure gradient highlights the fact that the model is only as robust as its underlying assumptions and will require additional refinement and revision as new experimental data on lens structure and function becomes available. In the following sections we discuss the assumptions and limitations of our current model with the view to highlight areas where further refinements of the model are required.

The solute flux in a fluid is governed by diffusion, electro-diffusion (if the solute is charged) and advection. Diffusion is the random walk of particles due to Brownian motion. Electro-diffusion is the flux of a charged particle due to the force applied by an electric field. Advection is the transport of a solute by a fluid that is moving. The physical origins of these transport processes and the derivations of the equations have been discussed previously in the literature [31, 32]. The implementation of these equations in a computational platform is very well explained elsewhere. In this study, we implemented the driven equations and solved for a set of converged-upon 3D fields. This method produced the calculated standing ionic concentration gradients in the lens predicted by our model [Figure 4.These gradients in turn gave rise to a trans-membrane electric-potential gradient field [Figure 5, based on the Hodgkin–Huxley model. Consider a cell membrane that is not equally permeable to all the present ionic species on either side of it. Here, the permeable ions will tend to move down their concentration gradient taking their electrical charge with them as they go. Therefore, an electrical potential will be generated, which will drive them in the reverse direction. An electrochemical equilibrium will be reached when the diffusive force equals the electromotive force.

Our model has calculated a hydrostatic pressure gradient [Figure 6, the existence of which has recently been confirmed experimentally [29]. While the orientations of the calculated and measured pressure gradients are similar, they differ in magnitude (48kPa compared to 10kPa). It has been proposed that this pressure gradient is generated by the restricted flow of water from the centre to the periphery of the lens by a pathway mediated by gap junction channels; since genetic manipulations to increase or decrease gap junction numbers produces inverse changes in the pressure gradient. This illustrates that structural components of the lens can influence the magnitude of the pressure gradient. This suggests that the difference between calculated and measured pressure fields may reflect the absence of a structural feature not currently captured in our model. In this regard we have recently identified a zone in the inner cortex of the lens that exhibits a reduction in the penetration of solutes and water [19, 22] that could influence the magnitude of the calculated pressure gradient. Future updates of the model will include such newly discovered structural elements, allowing their effect on the calculated pressure fields to be assessed.

The electrochemical fields, combined with the hydrostatic pressure gradients [Figure 6 in our model generate the circulating ionic currents throughout the lens [Figure 8 & Figure 9. The ionic fluxes are accompanied by water flows in the microcirculation system. The water fluxes are generally described by the Navier–Stokes equations, which are derived from the conservation of mass, momentum, and energy principals. Also, the fluid flow in the lens can be described as slow or low-Reynolds number flow. It is also reasonable to assume that the fluid flow in a normal lens is near or at steady-state at all times [18]. These simplifications reduce the Navier–Stokes equations to the Stokes equations. The derivation and implementation of these equations into a computational framework is discussed elsewhere [16, 18]. In our model, solving this set of equations results in the calculation of the water flow velocity fields [Figure 10. Although the current model appears to accurately mimic the physiological homeostasis of the lens [Table 6 & Table 7, there are still some aspects that need future improvement.

Our model presently has been implemented with a line suture structure which runs from the anterior to the posterior of the lens, through its core [Figure 3A&B]. However the mouse has a Y-shaped suture that rotates 180 °C from the anterior to the posterior pole [42]. Hence, we are planning to improve the current model with an anatomically accurate asymmetrical 3D structure of the sutures. In our current model, the 3D extracellular solute diffusion coefficients are constant throughout. However, we [21, 53] and others [54–56] have recently identified a barrier to the movement of solutes in the lens, using variety of techniques. We strongly believe that the existence of this barrier is very important in shaping the fluid dynamics of the lens and in general its microcirculation. Hence, we will implement the imaged barrier as a part of continuous improvement of our computational model.

Another important structural feature of the ocular lens is the presence of a Gradient of Refractive Index (GRIN), which acts to correct for inherent spherical aberration to improve the optical properties of the lens [57]. This GRIN profile of the lens is directly dependant on the local water/protein concentration makeup of the lens [58, 59]. Recently, we have experimentally showed that the water/protein gradient of the lens is actively upheld by the microcirculation system [19]. Hence, another step in the improvement of the current model is to add equations to estimate the concentration gradient of water in the lens. Using that calculated field, we will be able to estimate the 3D GRIN maps of the lens. We would then use these 3D GRIN maps and optical ray-tracing software to produce a model that links lens physiology to the optical properties of the lens.

Our first generation 3D finite element model of lens structure and function describes ion and fluid dynamics in the mouse lens. We chose to model the mouse lens as ion and fluid dynamics have been extensively studied in this species [3, 4, 15, 16]. We also believe the model is an essential first step towards creating a comprehensive model of the human lens. Any model of the human lens model would need to include its more complex structural features and would need to be created so that its dimensions could be altered to study the effects of lens growth and ageing on the circulation system [60]. This future model would enable us to study the changes in lens physiology thought to underlie the initiation of age related nuclear cataract.

## Conclusion

During this project, a 3D model of the flux movements inside the mouse ocular lens was designed and executed using our high performance computer (HPC). Reviewing the results of the current model, it appears that solute fluxes, accompanied by water, enter the lens via the extracellular space all around it but with larger magnitudes around the polar regions. Among these inwardly extracellular solute fluxes, the Na^{+} fluxes were seemed to be dominant followed closely by Cl^{-} fluxes. Conversely, the solute effluxes appear to be via the intracellular space and seemed to be more pronounced around the equatorial region of the lens. The K^{+} fluxes were found to be the primary intracellular fluxes, caused mainly by exterior Na^{+}/K^{+} ATPase pumps. The net effect of these influx and effluxes were thought to be best explained by the calculated net current densities. The pattern of these net current densities at the surface of the lens was similar to previous experimental findings
[5, 6, 52]. Same fields modelled inside the lens found to be in agreement with the microcirculation theory of the lens
[3, 9, 16, 61].

This study brings together all the available experimental and theoretical data on the fluid dynamics of the ocular lens in order to create a comprehensive 3D model of this tissue. Previous studies have investigated the links between the steady state fluxes in the lens and its physiological homeostasis [62–65]. Using our computational model, we would be able to study the connections between the biodynamic natures of these perturbations and their functional consequences.

## References

Bassnett S:

**Lens organelle degradation.***Exp Eye Res*2002,**74:**1–6. 10.1006/exer.2001.1111Mathias RT:

**Steady-state voltages, ion fluxes, and volume regulation in syncytial tissues.***Biophys J*1985,**48:**435–448. 10.1016/S0006-3495(85)83799-1Mathias RT, Rae JL, Baldo GJ:

**Physiological properties of the normal lens.***Physiol Rev*1997,**77:**21–50.Donaldson P, Kistler J, Mathias RT:

**Molecular solutions to mammalian lens transparency.***News Physiol Sci*2001,**16:**118–123.Robinson KR, Patterson JW:

**Localization of steady currents in the lens.***Curr Eye Res*1982,**2:**843–7. 10.3109/02713688209020020Patterson JW:

**Characterization of the equatorial current of the lens.***Ophthalmic Res*1988,**20:**139–42. 10.1159/000266570Tamiya S, Dean WL, Paterson CA, Delamere NA:

**Regional distribution of Na, K-ATPase activity in porcine lens epithelium.***Invest Ophthalmol Vis Sci*2003,**44:**4395–9. 10.1167/iovs.03-0287Mathias RT, Rae JL, Eisenberg RS:

**Electrical properties of structural components of the crystalline lens.***Biophys J*1979,**25:**181–201. 10.1016/S0006-3495(79)85284-4Mathias RT, Rae JL:

**Steady state voltages in the frog lens.***Curr Eye Res*1985,**4:**421–430. 10.3109/02713688509025156Mathias RT, Rae JL, Ebihara L, McCarthy RT:

**The localization of transport properties in the frog lens.***Biophys J*1985,**48:**423–434. 10.1016/S0006-3495(85)83798-XMathias RT, Rae JL, Eisenberg RS:

**The lens as a nonuniform spherical syncytium.***Biophys J*1981,**34:**61–83. 10.1016/S0006-3495(81)84837-0Rae JL, Kuszak JR:

**The electrical coupling of epithelium and fibers in the frog lens.***Exp Eye Res*1983,**36:**317–326. 10.1016/0014-4835(83)90114-8Zampighi GA:

**Distribution of connexin50 channels and hemichannels in lens fibers: a structural approach.***Cell Commun Adhes*2003,**10:**265–270.Candia OA, Zamudio AC:

**Regional distribution of the Na+ and K+ currents around the crystalline lens of rabbit.***Am J Physiol Cell Physiol*2002,**282:**C252.Mathias RT, White TW, Gong X:

**Lens gap junctions in growth, differentiation, and homeostasis.***Physiol Rev*2010,**90:**179–206. 10.1152/physrev.00034.2009Mathias RT, Rae JL:

**Transport properties of the lens.***Am J Physiol Cell Physiol*1985,**249:**181–190.Mathias RT, Wang H:

**Local osmosis and isotonic transport.***J Membr Biol*2005,**208:**39–53. 10.1007/s00232-005-0817-9Malcolm DTK:

*A Computational Model of the Ocular Lens [Internet]*. ResearchSpace@ Auckland, Auckland, New Zealand; 2006. [cited 2012 Sep 11]. Available from: https://researchspace.auckland.ac.nz/handle/2292/754Vaghefi E, Pontre BP, Jacobs MD, Donaldson PJ:

**Visualizing ocular lens fluid dynamics using MRI: manipulation of steady state water content and water fluxes.***Am J Physiol Regul Integr Comp Physiol*2011,**301:**R335-R342. 10.1152/ajpregu.00173.2011Vaghefi E, Jacobs MD:

**Uptake and distribution of gadolinium in the ocular lens.***Conf Proc IEEE Eng Med Biol Soc*2008,**2008:**843–846.Vaghefi E, Walker K, Pontre BP, Jacobs MD, Donaldson PJ:

**Magnetic resonance and confocal imaging of solute penetration into the lens reveals a zone of restricted extracellular space diffusion.***Am J Physiol Regul Integr Comp Physiol*2012,**302**(11):R1250–1259. 10.1152/ajpregu.00611.2011Vaghefi E, Pontre B, Donaldson PJ, Hunter PJ, Jacobs MD:

**Visualization of transverse diffusion paths across fiber cells of the ocular lens by small animal MRI.***Physiol Meas*2009,**30:**1061. 10.1088/0967-3334/30/10/007Vaghefi E:

*Computational modeling and magnetic resonance imaging of microcirculation in the ocular lens*. (Bioengineering)--University of Auckland, Auckland, New Zealand; 2010.Vaghefi E, Osman NAA, Abas WABW, Wahab AKA, Ting H-N:

**Ocular Lens Microcirculation Model, A Web-Based Bioengineering Educational Tool.**In*5th Kuala Lumpur International Conference on Biomedical Engineering 2011*. Springer Berlin Heidelberg, Volume 35. Berlin, Heidelberg; 2011:25–28.Vaghefi E, Hunter PJ, Jacobs MD:

**3D Finite Element Modeling of Avascular Circulation in the Ocular Lens.**In*4th Kuala Lumpur International Conference on Biomedical Engineering 2008*. Edited by: Abu Osman NA, Ibrahim F, Wan Abas WAB, Abdul Rahman HS, Ting H-N. Springer Berlin Heidelberg, Berlin, Heidelberg; 2008:469–472.Duncan G:

**Permeability of amphibian lens membranes to water.***Exp Eye Res*1970,**9:**188–197. 10.1016/S0014-4835(70)80075-6Delamere NA, Duncan G:

**A comparison of ion concentrations, potentials and conductances of amphibian, bovine and cephalopod lenses.***J Physiol*1977,**272:**167–186.Guerschanik SN, Reinach PS, Candia OA:

**Chloride compartments of the frog lens and chloride permeabilities of its isolated surfaces.***Invest Ophthalmol Vis Sci*1977,**16:**512–520.Gao J, Sun X, Moore LC, White TW, Brink PR, Mathias RT:

**Lens intracellular hydrostatic pressure is generated by the circulation of sodium and modulated by gap junction coupling.***J Gen Physiol*2011,**137:**507–520. 10.1085/jgp.201010538Benedek GB, Villars FMH

**2nd ed. Springer.***Physics with Illustrative Examples from Medicine and Biology: Statistical Physics*2000.Benedek GB, Villars FMH

**2nd ed. Springer.***Physics With Illustrative Examples from Medicine and Biology: Electricity and Magnetism*2000.Fischbarg J, Diecke FP, Kuang K, Yu B, Kang F, Iserovich P, Li Y, Rosskothen H, Koniarek JP:

**Transport of fluid by lens epithelium.***Am J Physiol*1999,**276:**C548–557.Candia OA, Gerometta R:

**Fluid movement across the surface of the isolated bovine lens.***Invest Ophthalmol Vis Sci*2003,**44:**3455.Candia OA:

**Electrolyte and fluid transport across corneal, conjunctival and lens epithelia.***Exp Eye Res*2004,**78:**527–535. 10.1016/j.exer.2003.08.015Varadaraj K, Kumari S, Shiels A, Mathias RT:

**Regulation of aquaporin water permeability in the lens.***Invest Ophthalmol Vis Sci*2005,**46:**1393–1402. 10.1167/iovs.04-1217McLaughlin S, Mathias RT:

**Electro-osmosis and the reabsorption of fluid in renal proximal tubules.***J Gen Physiol*1985,**85:**699–728. 10.1085/jgp.85.5.699Benedek GB, Villars F:

*Physics, with Illustrative Examples from Medicine and Biology: Mechanics*. 2000.Cannell MB, Jacobs MD, Donaldson PJ, Soeller C:

**Probing microscopic diffusion by 2-photon flash photolysis: measurement of isotropic and anisotropic diffusion in lens fiber cells.***Microsc Res Tech*2004,**63:**50–57. 10.1002/jemt.10422Gao J, Sun X, Yatsula V, Wymore RS, Mathias RT:

**Isoform-specific function and distribution of Na/K pumps in the frog lens epithelium.***J Membr Biol*2000,**178:**89–101. 10.1007/s002320010017Jacobs MD, Soeller C, Sisley AMG, Cannell MB, Donaldson PJ:

**Gap junction processing and redistribution revealed by quantitative optical measurements of connexin46 epitopes in the lens.***Invest Ophthalmol Vis Sci*2004,**45:**191–9. 10.1167/iovs.03-0148Parmelee JT:

**Measurement of steady currents around the frog lens.***Exp Eye Res*1986,**42:**433–41. 10.1016/0014-4835(86)90003-5Kuszak JR, Mazurkiewicz M, Zoltoski R:

**Computer modeling of secondary fiber development and growth: I.***Nonprimate lenses. Molecular Vision*2006,**12:**271–282.Fischer G, Tilg B, Modre R, Huiskamp G, Fetzer J, Rucker W, Wach P:

**A bidomain model based BEM-FEM coupling formulation for anisotropic cardiac tissue.***Ann Biomed Eng*2000,**28:**1229–1243.Duncan G, Jacob TJ:

**Influence of external calcium and glucose on internal total and ionized calcium in the rat lens.***J Physiol*1984,**357:**485–493.Zampighi GA, Eskandari S, Kreman M:

**Epithelial organization of the mammalian lens.***Exp Eye Res*2000,**71:**415–435. 10.1006/exer.2000.0895Blankenship T, Bradshaw L, Shibata B, FitzGerald P:

**Structural specializations emerging late in mouse lens fiber cell differentiation.***Invest Ophthalmol Vis Sci*2007,**48:**3269. 10.1167/iovs.07-0109Duncan G:

**Movement of sodium and chloride across amphibian lens membranes.***Exp Eye Res*1970,**10:**117. 10.1016/S0014-4835(70)80017-3Hunter PJ, Borg TK:

**Integration from proteins to organs: the Physiome Project.***Nat Rev Mol Cell Biol*2003,**4:**237–243. 10.1038/nrm1054Rubin J, Wünsche BC, Cameron L, Stevens C:

**Animation and modelling of cardiac performance for patient monitoring.***In Proceedings of IVCNZ*2005,**5:**476–481.Wang H, Gao J, Sun X, Martinez-Wittinghan FJ, Li L, Varadaraj K, Farrell M, Reddy VN, White TW, Mathias RT:

**The effects of GPX-1 knockout on membrane transport and intracellular homeostasis in the lens.***J Membr Biol*2009,**227:**25–37. 10.1007/s00232-008-9141-5Delamere NA, Paterson CA:

**The influence of calcium-free solutions upon permeability characteristics of the rabbit lens.***Exp Eye Res*1979,**28:**45–53. 10.1016/0014-4835(79)90104-0Parmelee JT, Robinson KR, Patterson JW:

**Effects of calcium on the steady outward currents at the equator of the rat lens.***Invest Ophthalmol Vis Sci*1985,**26:**1343–8.Grey AC, Jacobs MD, Gonen T, Kistler J, Donaldson PJ:

**Insertion of MP20 into lens fibre cell plasma membranes correlates with the formation of an extracellular diffusion barrier.***Exp Eye Res*2003,**77:**567–574. 10.1016/S0014-4835(03)00192-1Moffat BA, Pope JM:

**Anisotropic water transport in the human eye lens studied by diffusion tensor NMR micro-imaging.***Exp Eye Res*2002,**74:**677–87. 10.1006/exer.2001.1164Truscott RJW:

**Age-related nuclear cataract-oxidation is the key.***Exp Eye Res*2005,**80:**709–25. 10.1016/j.exer.2004.12.007Moffat BA, Landman KA, Truscott RJ, Sweeney MH, Pope JM:

**Age-related changes in the kinetics of water transport in normal human lenses.***Exp Eye Res*1999,**69:**663–9. 10.1006/exer.1999.0747Smith G, Cox MJ, Calver R, Garner LF:

**The spherical aberration of the crystalline lens of the human eye.***Vision Res*2001,**41:**235–243. 10.1016/S0042-6989(00)00206-6Pierscionek BK, Chan DY:

**Refractive index gradient of human lenses.***Optom Vis Sci*1989,**66:**822–829. 10.1097/00006324-198912000-00004Ioseliani OR, PIERSCIONEK B

**illustrated. Nova Publishers.***Species variations in the refractive index of the eye lens. Focus on Eye Research*2006, 238.Dubbelman M, Van der Heijde G:

**The shape of the aging human lens: curvature, equivalent refractive index and the lens paradox.***Vision Res*2001,**41:**1867–1877. 10.1016/S0042-6989(01)00057-8Mathias RT, Kistler J, Donaldson P:

**The lens circulation.***J Membr Biol*2007,**216:**1–16. 10.1007/s00232-007-9019-yKobatashi S, Roy D, Spector A:

**Sodium/potassium ATPase in normal and cataractous human lenses.***Curr Eye Res*1982,**2:**327–334. 10.3109/02713688209000778Spector A:

**Oxidative stress-induced cataract: mechanism of action.***FASEB J*1995,**9:**1173.White TW, Goodenough DA, Paul DL:

**Targeted ablation of connexin50 in mice results in microphthalmia and zonular pulverulent cataracts.***J Cell Biol*1998,**143:**815–25. 10.1083/jcb.143.3.815Buzhynskyy N, Girmens JF, Faigle W, Scheuring S:

**Human cataract lens membrane at subnanometer resolution.***J Mol Biol*2007,**374:**162–169. 10.1016/j.jmb.2007.09.022

## Acknowledgments

We wish to thank Prof Richard T Mathias for his insightful inputs in the development of the current model.

## Author information

### Authors and Affiliations

### Corresponding author

## Additional information

### Competing interests

The authors declare that they have no competing interests.

### Authors' contributions

**EV** created the computational model and drafted this manuscript. **DM** checked the integrity of the equations and numerical presentations. **MJ** checked the numerical outcomes of the model and edited the paper. **PD** conceived the manuscript; final edited it and approved the final version. All authors read and approved the final manuscript.

## Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

## Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## About this article

### Cite this article

Vaghefi, E., Malcolm, D.T., Jacobs, M.D. *et al.* Development of a 3D finite element model of lens microcirculation.
*BioMed Eng OnLine* **11**, 69 (2012). https://doi.org/10.1186/1475-925X-11-69

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/1475-925X-11-69

### Keywords

- Computational modelling
- Ocular lens
- Microcirculation
- Finite element
- Physiological optics