Mathematical and Numerical Methods for Modeling Drug Delivery to the Posterior Segment of the Eye

The high incidence of diseases that affect the posterior segment of the eye (PSE) here intended as composed by sclera, choroid and retina-prompts for establishing effective and well tolerated therapies. Topical application (instillation of drops) and systemic assumption remain the most widespread drug administration routes. However, the drug achieved levels are not therapeutically sufficient therapeutic, since in the first case the drug is mainly washed away by different pathways (aqueous humor, systemic adsorption, tears) and in the second reaches the PSE target in a minimal fraction. Intravitreal injections and biodegradable episcleral implants have emerged in the last decades as alternative, more effective, administration routes. Whilst these techniques already offer significant improvements, much space is still open to research. The difficulty in delivering drugs to the PSE via the intravitreal (IV) or episcleral (EP) route resides in the several physical and dynamical barriers-including blood retinal barrier, clearance from choriocapillaries and lymphatics-which hinder the passage of the drug molecules. In combination, and in support of clinical experiments, mathematical and computational methods can be used to simulate drug concentration levels in the tissues upon IV or EP administration. In this respect, this short review gives an update of the most recent (from year 2000 to present time) developments on mathematical models for drug delivery to the PSE. We specifically focus our attention on physiological modeling works that include spatial dependency. Our review work is organized in short sections accompanied by detailed tables. We discuss descriptions of the PSE morphology, considering “anatomically accurate” as well as reduced models and we analyze the biophysical phenomena included in the examined models. We present the numerical techniques adopted to solve the resulting systems of partial differential equations and we deal with the delicate issue of parameter choice and validation of the results. Eventually, we examine the main results and scientific achievements of the considered models. Our conclusions point out the absence of a “standard mathematical model” and highlight a significant scattering of the results obtained from the different authors.


Introduction
The unique anatomy and physiology of the eye makes drug delivery to the PSE -here intended as composed by sclera, choroid and retina -a highly challenging task [1]. In the past, drugs targeting the PSE were administered for the most part via topical or systemic routes, but severe limitations and inefficiency sources were associated with these therapeutic approaches. Intravitreal injections and (juxta-) episcleral implants have more recently emerged important tools of analysis. Previous reviews specifically focusing on mathematical modeling of drug delivery to the PSE can be found in [4,5]. We refer to [2,[6][7][8][9] for upto-date, comprehensive, reviews on physiological and clinical aspects of this topic.
The present short review gives an update of the most recent (from year 2000 on) developments on mathematical modeling of drug delivery to the PSE. We specifically focus our attention on physiological modeling works that include spatial dependency in the interested domains. This choice stems from the above observations about the inhomogeneous distribution of drug in the complex PSE structures, which definitely prompts for including spatial dependency in the models. We do not include mathematical models tailored to study drug distribution in the sole vitreous, for which a specific and relatively vast literature exists, see [10,11]. In our analysis, we have identified 11 papers which as alternative routes to achieve sustained therapeutic drug levels [2]. The penetration rate of a certain drug through these latter routes is influenced by its solubility, hydro/lipophilicity, charge, degree of ionization and molecular size [3]. Furthermore, drug concentration in the PSE tissues is not homogeneous but distributes in a complex way according to several concurrent biophysical phenomena as convection, diffusion, losses through anterior and posterior pathways, clearance mechanisms through blood and lymphatic vessels and degradation processes (see also Table 1 of [4]. It has also been observed [4] that only a weak correlation exists between the actual local drug concentration in the PSE, which determines the pharmacological response, and the mean concentration values that experiments can monitor in the more accessible aqueous chamber and vitreous body under in vivo conditions. In this perspective, mathematical and computational models can offer main results and scientific achievements of the models. Eventually, in Sect. 6 we draw the conclusions.

Geometrical Modeling
Morphometric information on the PSE geometry has been traditionally obtained from dissection of specimina and, more recently, from digital imaging, both in animals and humans. In establishing a geometrical model for the PSE, important issues must be faced: i) which are the relevant anatomical structures to consider as mathematical domains of the model? ii) What are anatomically correct measures of the included volumes? (for example, the volume of the vitreous fulfill the above described criteria, and whose content is synthetically reviewed in table 1, table 2 and table 3. Each  table is focused on a particular aspect of mathematical modeling of drug delivery to the PSE: representation of the PSE geometry, modeling of the relevant biophysical phenomena and discretization techniques adopted to obtain the numerical solution of the model on a computer. We organize the review as follows. In Sect. 2, we analyze and compare the geometrical modelling choices; in Sect. 3, we present the biophysical phenomena included in the models; in Sect. 4, we deal with the numerical techniques adopted to solve the models, and in Sect. 5 we discuss the the biophysical phenomena included in the considered mathematical models. All the models include diffusion mechanisms for the drug and the majority also includes porous/fluid filtration, modeled for the most part by Darcy equations, but also by Navier-Stokes equations [15] and by Darcy-Brinkman equations [14]. Simulation results show that the filtration field has almost a negligible effect for low-weight drugs [17,18], whilst it becomes significant for high-weight drugs. On the contrary, the contribution of RPE pumping always appear important, especially when the retina is the therapeutic target. The RPE pumping effect, included in a very limited selection of papers, is modeled as a fictitious advective field or through different inward/outward permeabilities of the RPE membrane [13,19,20]. The inclusion of the significant clearance pathways for a certain drug administration mode is rather varied. Choriocapillaries are often treated as a perfect sink [16,21] or they are modeled via porous jump equations [13]. Retinal blood vessels are only considered in [18] and, via a mixture equation approach, in [20]. Generally, the fluid filtration is assumed to be steady state (except for) [15] and not influenced by the drug transport.

Numerical Approaches and Validation Techniques
In table 3, we provide a comparative analysis between the molecules considered in each model, the mathematical description of the drug source, the numerical approaches and parameter estimates used in the considered models for drug delivery to the PSE. For the most part, lowweight molecules are considered, with a predominance of the marker fluorescein on which the most part of experimental data are based. Only a few papers address different molecules with therapeutic significance, as the anti-VEGF macromolecular drug ranibizumab (IgG1 Fab) in [17]. The drug source is represented for IV injection as a bolus of prescribed concentration of spherical shape [15,21] or time evolving shape [12] and for EP plug as a prescribed concentration [13] or prescribed time dependent concentration [17] at the scleral boundary. In Ferreira, et al. [22], the drug implant is inserted in the vitreous.
The inclusion in the model of the spatial dependency leads to large systems of partial differential equations, whose solution with analytical techniques is unfeasible, if not impossible. Thus, the use of numerical methods is in order. Spatial discretization can be based on finite elements (FEMs) [13,19] or finite volume techniques [15], generally implemented in commercial packages or on finite differences [18]. Different discretization schemes are adopted for time marching, including implicit/semi-implicit finite difference schemes [13,14] and adaptive Runge-Kutta schemes [20]. Three-dimensional problems give rise to large size systems, with a number of degrees of freedom ranging from about 200,000 nodes [13] to 321,000 chamber considered in Jooybar, et al. [12] is equal to 4.4 ml , where as in the original model of [13], the volume of the vitreous is about 2.2 ml); iii) what are the significant clearance pathways (for example, blood vessels, periocular structures or lymphatic vessels) for a certain drug administration mode? iv) how to represent un-resolved structures via boundary conditions on the external surfaces? How to deal with internal interface/ coupling conditions? In addition, one should consider the fact that the PSE has unique geometrical features in each subject. In this respect, the common procedure consists in extensive literature searches and data comparisons so to obtain models based on a representative, average, geometry for a certain species (human, rabbit, mouse, monkey.
In table 1, we provide a comparative analysis between the geometrical models of the PSE adopted in the considered mathematical works (being for the totality 3D or 1D models, except for the 2D model of [14]. The most part of the papers considers geometrical parameters corresponding to the human eye, possibly partially derived from other species (rabbit eye, see also the comparative analysis carried out in [15]. Moreover, for the most part, these recent models consider, at least, vitreous and retina as separated domains. Choroid and sclera are sometimes lumped together in a unique domain, as in [16]. Anterior structures are described with different levels of precision, ranging from the "anatomically accurate" representation of [15] to the representation of mimicking boundary conditions of [17]. Interface conditions between adjacent layers l and m are, in their most general form, expressed by the permeation relation L P C C n ∂ = − ∂ D l being the diffusion coefficient in layer l with outward unit normal vector n l , C l and C m the drug concentrations in layers l and m, respectively, L lm the porous permeability of the membrane between layers l and m and P lm the partition coefficient keeping into account the affinity of the drug in the two layers. As a general remark, it is worth noticing that, while in [4] it is emphasized that fine geometric details are not needed to capture the salient features of drug distribution and only the proportions of different domains of the eye are considered to be important, more accurate analyses may need the inclusion of structures responsible of specific pathways (as for example the periocular space in [18] or the gap of Petit in [15].

Biophysical Phenomena Modeling
In table 2, we provide a comparative analysis between from analytical relations/correlations between different layers [18]. It should be observed that, whilst the actual value of some of these parameters will not probably affect too much the results when dealing with low weight molecules, when considering larger/heavier molecules the parameter values can have a more significant impact [21].

Results
In table 4, we present the main results obtained in the considered mathematical models. Simulations have been performed to study concentration profiles in different conditions. Sensitivity studies have been carried out to address the significance of: i) certain biophysical mechanisms (active transport in the retina) [13,20], clearance rates [15]; ii) technical details in the drug administration (needle angle) [12], implant position [14]; different boundary conditions (vitreous outflow) [21]; implications of interspecies variability [15]; pathological conditions (inflammation of the retina), posterior vitreous detachment [22].
nodes [17] and may require the use of parallel computing (as in the simulations carried out on a quad processor in [15]. Special care must be taken when simulating diffusion/convection/reaction equations for drug delivery. Convection dominates when considering high weight/ scarcely diffusing molecules or modeling RPE pumping actions as a fictitious velocity field. Reaction, representing metabolic terms/degradation effects, may also become dominant and induce instabilities and/or unphysical negative concentration fields [4]. Best fit between the simulated drug levels and experimental/in vivo data have been performed to assess the model validity and fit critical parameters (using for the most part fluorescein data as those from Palestine & Brubaker) [23] in Balachandran & Barocas [13] or those from [24] in Park, et al. [21]; Gd-DTPA data were used in [16]. Parameters estimated from data fitting of fluorescein have been often used also for different drugs/delivery routes. Diffusivity, permeability and RPE pumping parameters have been computed from data fitting [13,16,18], from analysis of TEM images [14] or Table 4: Validation method and main results for 3D physiological models of drug delivery to the PSE.

Reference
Validation method Main results [16] Gd-DTPA data from MRI in vivo and ex vivo estimate of lumped posterior membrane (6.0 × 10 -8 cm 2 /s) and vitreous (2.8 × 10 -6 cm 2 /s) diffusivity coefficient [21] fluorescein data from in vivo IV injection comparison between IV injection and implant; influence of vitreous drug diffusivity as a surrogate for drug molecular weight; influence of changes of the proportion of fluid loss by vitreous and aqueous outflows [18] validation from GFP data in mouse after periocular injection estimation of resistances (RPE: 9 × 10 5 , periocular: 1.6 × 10 -1 and choroid: 3.5 × 10 -3 ), estimation of clearance rates from fitting of experimental data; sensitivity analysis on resistance, clearance rates, diffusivities, void volume for peak concentration, peak time and total integrated concentration [13] validation from experimental data for systemic administration study of the influence on drug levels of RPE active pumping, episcleral losses and choriocapillary clearance [15] comparison with experimental data in rabbit eye after intravitreal injection estimate of the hydraulic resistance of the trabecular meshwork by fitting procedure; comparison of solutions (fluid flow and drug) in different geometries; parametric study of clearance from the ASE; methodology to excerpt results for one species (human, rabbit, monkey) from different species [14] transient mean plasma concentration and contours of anecortave acetate for different blood velocities; mean concentration in the choroid as indicator of the drug bioavailability at the retina [17] validation from fluoroscein data (IV injection, systemic) effect of spatial inhomogeneities on episcleral drug delivery (to be compared with compartmental model by) [27], estimate of time duration of therapeutic drug levels for the igG1 Fab molecule [19] effect drug mobility in sclera due to diffusion and permeability and of the partition coefficient at retina/vitreous interface and at retina/choroid interface [12] validation from experimental data for systemic (intravenous) administration study of the effect of injection time and needle gauge/angle; study of the effect of implant position/type [20] comparison with results obtained from recent models [27,13,18] sensitivity analysis on diffusivity, clearance rates, permeability coefficient for drug concentration; sensitivity analysis on filtration velocity and fictitiousactive velocity at RPE; theoretical analysis on lower and upper bound of drug concentration in the retina [22] influence In general, experimental data are very sparse and significantly differ from each other. This lack of "golden standard experimental data" makes more difficult to dispose of well assessed data for model validation purposes 5. The mathematical representation itself of the phenomena (for example, representing RPE pumping effects by a fictitious convective field or by an altered permeability) impacts on the results. However, again, there is no "gold standard" for carrying out a comparison.

Conclusions and Perspectives
According to the World Health Organization, in 2010, more than 13% of the cases of blindness were the result of diseases of the PSE [25]. It is thus of high importance to dispose of efficient and well tolerated ways to deliver drugs to this region. Mathematical models can provide tools useful to enlarge the knowledge on this topic, supporting and completing clinical experiments. In general, much space is still open when it comes to address therapeutically relevant drugs, especially macromolecules (anti-angiogenic antibodies, oligonucleotides, growth factors, and trans-gene expression) [4] and prodrugs improving drug bio-availability) [6]. For these drugs, parameter estimation may need reviewing along with reconsideration of the significant biophysical With the purpose of comparison of the different models when computing a similar quantity, in figure 1 we show drug peak (and mean) concentration value and time-to-peak in the retina and choroid. We represent the data where available. The results are all obtained for lowweight molecules comparable to fluorescein. The drug is delivered via EP plug or EP injection. It is evident from the graphs a significant scattering of the results, of which several factors may be responsible. What appears more evident from our review work is that: 1. There is no general agreement on the detail level required for the geometrical description of the anatomical structures. Different choices are carried out, depending on the considered drug delivery route (IV injection, episcleral implant and even with the same approach with the precise location of the drug source), species (man, rabbit, monkey, rat) and on the intraspecies variability 2. The characteristics of the considered drug (hydrophilicity/lipophilicity, molecular weight and size, charge) can impact the choice of including or not in the model certain bio-physical phenomena (for example, the clearance routes and the permeability across membranes as the RPE) for a specific drug. Moreover they can also influence, by themselves, the geometrical detail needed in the description of a certain structure 3. The value of the parameters strongly influence the model results. This is clearly demonstrated by the results obtained by the sensitivity analyses performed in several  [19]. Values from [13] have been interpreted as percentage of the enforced boundary concentration value, in lack of indications. mechanisms. To cite one example, the convection field can become a determinant of the drug fate, and contrast/ enhance RPE pumping actions. A better understanding of the impact of aging/pathological conditions is necessary to develop more effective therapies so that drug concentration is maintained within the therapeutic window at the target site for the desired period of time (so far the effect of the liquefaction of the vitreous, breakdown of blood retinal barrier and elevated intraocular pressure have been considered in Ferreira, et al.) [22]. This aspect should also include the effect of surgical procedures, for example the increased drug clearance due to vitrectomy [8] or the modified fluid-dynamics due to the use of vitreous substitutes [26].