Development of a 3D finite element model of lens microcirculation

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.

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][9][10][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  (EC), that at the equator (EQ) divide to produce the elongating differentiating fibre cells (DF), that eventually lose their nuclei and cellular organelles to become mature fibre cells in the lens nucleus (MF). Fiber cells from adjacent hemispheres meet at the anterior (AP) and posterior (PP) poles to form the sutures. Arrows in the diagram represent the direction of ion and water fluxes. These fluxes have been directly measured outside the lens [5,6] but their position and direction inside the lens are to date purely theoretical [3]. B: An equatorial cross section of the lens showing a cellular view of ion and water movement in the lens [7]. Current and solutes are proposed to flow into the lens via the extracellular space, to cross fibre cell membranes, and to flow outward via an intracellular pathway mediated by gap junction channels.
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 Clconductances [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][16][17][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][20][21][22][23][24][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][27][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][16][17][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][27][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][33][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][16][17][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  Table 3] to account for the tortuosity of the extracellular cleft. The conductivities for the Na + , K + and Clchannels 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 electroosmosis 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].
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. ]. 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 Figure 2 Steps of the iteration solution of the model reviewed in the text. The first step solved the Stokes equations, using the initial conditions (C 0 , φ 0 ) and boundary conditions (C αo , φ o , p o ), to calculate the fluid velocity (u n+1 ) and pressure fields (p n+1 ). The second step used electro-neutrality to calculate new potential fields (φ n+1 ).
Step three calculated new solute fluxes (j n+1 ) using current concentration fields (C n ), new potential field (φ n+1 ) and new velocity field (u n+1 ). The forth step used current concentration fields (C n ) and new potential field (φ n+1 ) to calculate the new trans-membrane or surface solute sources (s n+1 ).
Step five involved using the newly calculated solute sources (s n+1 ) to calculate new concentration fields (C n+1 ). The model then checked for the convergence criterion automatically and updated the initial fields and repeated these steps if the criterion was not met. 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][27][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 (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 Clin 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 Clin the lens have also been measured using chemical analysis in a number of species [26][27][28]47]. To facilitate comparison between these previous measurements of total lens ion concentrations and the values Note that in the extracellular space, a diffusiveconductive gradient was created by the decreasing concentration of Na + from the periphery towards the centre. In the intracellular space, however, the concentration field for Na + was reversed. The small extraand intracellular Clconcentration field in the radial direction is close to its predicted steady state [3]. K + is transported into cells mostly close to the periphery of the ocular lens [16] and so its high concentration at the intracellular space of the surface and accumulation in the core of the extracellular space. The values on the colour-bars are in Molar units. 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 Clappeared 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].
(See figure on previous page.) Figure 5 Standing electrical fields modelled in the mouse lens. 3D quarter section views of the standing electrical fields calculated from the model for the extracellular (A) and intracellular (B) space plus the transmembrane potential field (C). Both the extracellular and intracellular electrical fields decrease in value from surface to centre, but the magnitude of the observed change in the intracellular potential is smaller due to the high level of gap junction coupling between inner-cells and surface cells. The transmembrane electrical potential was calculated as the difference between the intracellular and extracellular potentials. This field appeared to be decreasing in magnitude from the surface towards the core of the lens. The values on the colour-bars are in Volt units.  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].   (A, B), Cl -(C, D) and K + (E, F). Extracellular Na + fluxes seemed to be largest of the ionic fluxes with a maximum located around the polar regions of the model. The intracellular fluxes on the hand were outwardly and equatorially oriented. Extra and intracellular Clfluxes appeared to follow the same path into the lens with close magnitudes to Na + fluxes in the same compartment. The K + extracellular fluxes seemed to be negligible compared to other ions, while intracellular values were the largest and highly concentrated around the equatorial region of the model. The numbers on the colour-bar are in (mMol/cm 2 .s) units.

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 electricpotential 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][55][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 raytracing 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 Clfluxes. 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][63][64][65]. Using our computational model, we would be able to study the connections between the biodynamic natures of these perturbations and their functional consequences.