A shell model of eye growth and elasticity
L. S. Kimpton, B. J. Walker, C. L. Hall, B. Bintu, D. Crosby, H. M. Byrne, A. Goriely
aa r X i v : . [ phy s i c s . b i o - ph ] O c t Journal of Elasticity manuscript No. (will be inserted by the editor)
A shell model of eye growth and elasticity
L. S. Kimpton · B. J. Walker · C. L. Hall · B. Bintu · D. Crosby · H. M. Byrne · A. Goriely
Received: date / Accepted: date
Abstract
The eye grows during childhood to position the retina at the correct distance behind the lens toenable focused vision, a process called emmetropization. Animal studies have demonstrated that this growthprocess is dependent upon visual stimuli, while genetic and environmental factors that affect the likelihood ofdeveloping myopia have also been identified. The coupling between growth, remodeling and elastic responsein the eye is particularly challenging to understand. To analyse this coupling, we develop a simple model ofan eye growing under intraocular pressure in response to visual stimuli. Distinct to existing three-dimensionalfinite-element models of the eye, we treat the sclera as a thin axisymmetric hyperelastic shell which undergoeslocal growth in response to external stimulus. This simplified analytic model provides a tractable frameworkin which to evaluate various emmetropization hypotheses and understand different types of growth feedback,which we exemplify by demonstrating that local growth laws are sufficient to tune the global size and shapeof the eye for focused vision across a range of parameter values.
Keywords
Eye · Emmetropization · Myopia · Elastic shell · Morphoelasticity.
Mathematics Subject Classification (2020) · In health, the human eye grows during childhood in order to adopt the correct size and shape for focusedvision in a process called emmetropization . The goal of emmetropization is to position the retina at the correctaxial distance behind the lens for the optical power of the anterior eye. When this process fails, the individualis either myopic (short-sighted, with excessive axial length) or hyperopic (long-sighted, with insufficient axiallength). Whilst, for many, myopia is readily treatable with prescribed lenses, severe myopia (classified as arefractive error of − L. S. Kimpton · B. J. Walker · H. M. Byrne · A. GorielyMathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road,Oxford, OX2 6GG, UK.Tel.: +44 (0)1865 273525E-mail: [email protected]. L. HallDepartment of Engineering Mathematics, University of Bristol, Merchant Venturers Building, Woodland Road, Bristol,BS8 1UBD. CrosbyWave Optics Ltd, 41 Park Drive, Milton Park, Abingdon, OX14 4SR, UK.Eyejusters Ltd, Unit 6, Curtis Industrial Estate, North Hinksey Lane, Oxford, OX2 0LX, UK.B. BintuDepartment of Physics, Harvard University, Cambridge, Massachusetts 02138, USA. ncidence expected to continue to increase [4–6]. This is attributed to a number of factors, ranging from areduction in time spent outside to an increase in time spent reading [1, 7].There are treatments that aim to slow the progression of myopia in children, as discussed in a numberof comprehensive reviews [8–10], with the most effective being the use of atropine eye-drops, whilst othersinclude non-traditional corrective lenses. The theory underpinning the prescription of these lenses, whichunder-correct a myopic eye, assumes that myopic blur detected at the retina prevents or slows subsequentaxial growth, hence reducing myopia progression. Building upon this assumption, we suggest that a necessarystep towards understanding and improving the treatment of myopia is to first understand the process ofemmetropization. Thus, as the primary goal of this work, we will aim to develop a simple mechanical modelof eye growth that can be used to investigate hypotheses for emmetropization.The eye is an intricate organ, with many interconnected mechanical and optical components. We illustratesome of these features in figure 1, though in this work we will focus primarily on modeling the deformationand growth of the sclera, upon which the retina and underlying choroid sit, with minimal consideration ofthe anterior eye. The sclera contributes the majority of the mechanical stiffness of the eye, with its thicknessvarying from 0 . ScleraChoroidRetinaFoveaOptic nerveIrisCorneaLensLimbus
Fig. 1
A schematic diagram of a horizontal section through an eye, with many of the structural features labelled. In thiswork, we will restrict attention to explicit considerations of only the sclera, retina, and cornea, though we will also utilisethe optical properties of the anterior eye.
Indeed, something not enjoyed by complex numerical models is a suitability for analytical study andrapid exploration, a desirable trait of other models that has been successfully exploited in similar contexts toidentify plausible growth laws and instabilities [24–26]. These simpler models, often described as toy modelsor caricatures, utilise simplified geometries whilst retaining key properties of the full system, with the goalnot to obtain predictive models but to instead gain insight into the feedback mechanisms between growth,remodeling, geometry and elasticity. This general approach can be particularly valuable when consensus islacking about the underlying mechanisms, as is often the case in physiological systems. Thus, in the contextof emmetropization, there remains significant scope for a simplified modeling framework that captures themechanical and growth features of the eye, with such a framework enabling future evaluation and explorationof various hypotheses for the wide range of processes involved in ocular development.One such hypothesis concerns the drivers of scleral growth during emmetropization. Supported by variousobservations in animals [27–33], though perhaps not directly translatable to the human eye, it is suggested thatat least part of the growth stimulus derives from locally interpreted visual information on the retina, notablyin combination with a host of interacting genetic and environmental factors. Indeed, multiple works reportthat the axial length of the grown eye in guinea pigs and chicks is dependent on the wavelength of the incidentlight [29–31], suggesting a reliance of growth on the colour of observed light, which appears consistent withthe aforementioned incidence of myopia increasing in line with reduced time spent outdoors. Here, utilising2he optics of the anterior eye and posing a simple model for growth stimulus that phenomenologically capturesthis behaviour, we will consider this example hypothesis and test its ability to produce qualitatively realisticmorphologies.We will proceed by formulating a simple but versatile model of emmetropization, suitable for the explo-ration and evaluation of models of growth and mechanical features of the developing eye. In particular, wewill focus on the later stages of emmetropization, in which the optical properties of the anterior eye can beconsidered to be fixed, though this assumption may be relaxed in future study. In response to visual stimuli,which we will later define as the ability to correctly focus blue light at the retina, growth of the scleral shellwill act to relieve the local stress resulting from intraocular pressure, which we explore as a possible driver foremmetropization. We will also touch upon the effects of including sophisticated mechanical properties, such asfibre reinforcement due to collagen present in the sclera, as well as considering the effects of scleral thicknessso as to further showcase the tractability and flexibility of this model approach.
The sclera is modeled as a thin, morphoelastic shell that resists both bending and stretching and is inflatedby the intraocular pressure. Requiring the shell to remain axisymmetric under loading and growth permitsa choice of coordinates that allows the model to be formulated solely in terms of principal stretches, i.e. ourrepresentation of the deformation gradient is diagonal. To describe the growth and the elastic response ofthe sclera, we utilise the multiplicative decomposition approach formulated by [34], commonly referred to as morphoelasticity . In this formulation, the total deformation of the sclera in each principal direction is writtenas the product of an elastic stretch and a ‘growth stretch’. This has been applied in a wide range of biologicalapplications, as reviewed by [35] and [36] and described recently by [37].In our model, we suppose that the sclera deforms almost instantaneously in response to changes in in-traocular pressure, so that the elastic response is very fast compared to growth, which occurs on the timescaleof years. Consequently, we assume mechanical equilibrium of the sclera at each instant and determine theelastic stretches that characterise the deformation from the unloaded, grown, reference configuration to thecurrent, pressurised configuration. The current position of the sclera is used to estimate the degree of blurexperienced locally at the retina, assuming smooth anterior attachment to the cornea. This blur defines thegrowth stretches that are used to update or ‘grow’ the reference configuration at the next timestep. Themodel formulation is detailed below and closely follows that used previously to describe fungal growth andcell blebbing [38, 39].2.1 GeometryWe model the eye with a simplified geometry, with the centre of the pupil and the fovea both lying on theanterior-posterior axis, about which the sclera is axisymmetric. The physiological eye is not axisymmetric andthe fovea sits slightly temporal to the posterior pole, but axisymmetry is a reasonable simplification. Theposition of the sclera is defined by rotating the curve C , which represents the centreline of the sclera aboutthe z -axis (see figure 2), with the retina lying on the interior surface of the sclera. The curve is parameterised zs ( σ, t ) r ( σ, t ) C e s n θ ( σ, t )( z ⋆ , r ⋆ ) Fig. 2
Schematic diagram of the coordinate system in a plane of constant φ for the deformed and grown scleral shell. Theleftmost point is ( z ⋆ , r ⋆ ), the location of smooth attachment of the sclera to the cornea. by a material parameter σ ( Σ, t ), which is the arclength in the grown but unloaded configuration measuredfrom where C meets the z -axis at the rear of the eye. In turn, Σ ∈ [0 , L ] is an arclength parameter in theinitial, unloaded configuration, with smooth attachment to the cornea occurring at Σ = L . Returning to the3eformed configuration, s ( σ, t ) is the arclength distance from the z -axis, r ( σ, t ) is the radial distance to the z -axis and θ ( σ, t ) is the angle that the normal to C makes with the z -axis, as shown in figure 2, so that ∂r∂s = cos θ , (1a) ∂z∂s = − sin θ . (1b)We associate a unit outward normal e n and two unit tangent vectors e s and e φ with each point on the shell,which point respectively in the direction of increasing s and increasing φ , where φ is the angle that C is rotatedabout the z -axis. These tangent vectors e s and e φ are the principal directions associated with the principalcurvatures, κ s and κ φ , where κ s = ∂θ∂s , (2a) κ φ = sin θr . (2b)Finally, we denote the unloaded radius by ρ ( σ, t ), analogous to r ( σ, t ) in the absence of elastic deformation,so that the elastic stretches in the e s , e φ , and e n directions are α s ( σ, t ) = ∂s∂σ , (3a) α φ ( σ, t ) = rρ , (3b) α n ( σ, t ) = hH , (3c)respectively, where h is the deformed scleral thickness and H is the unloaded scleral thickness.2.2 Force and moment balancesInstantaneous mechanical equilibrium equations are obtained by balancing the forces and torques acting on asmall patch of the thin sclera as in [39]. The sclera deforms in response to ∆P , the difference in intraocular andambient pressure. We assume that both the intraocular and ambient pressure remain constant, so that ∆P isset as constant throughout the growth and remodeling process. The stress resultants t φ and t s act along theshell in the e φ and e s directions, respectively, and both have units of force per unit length. The shear stressresultant, which acts on surfaces with normal e s in the direction e n , is denoted by q s and has units of forceper unit length. We note that our assumption of axisymmetry ensures that no such shear acts on surfaces withnormal e φ . The bending moments about the e φ and e s directions are written as m s and m φ , respectively, andboth have units of force. Neglecting inertial effects, the momentum balance supplies ∂∂s ( rq s ) = r ( ∆P − κ φ t φ − κ s t s ) , (4a) ∂∂s ( rt s ) = κ s rq s + t φ cos θ , (4b) ∂∂s ( rm s ) = m φ cos θ − rq s . (4c)A discussion of constitutive assumptions, used to define t s , t φ , m s and m φ and close these mechanical equations,is postponed to section 2.4, along with appropriate boundary conditions as derived in section 2.6.2.3 Morphoelastic stretches and growthAs outlined above, we assume a multiplicative decomposition of the deformation, with the total deformationfrom the initial, unloaded configuration to the grown, loaded configuration defined to be the product of an4lastic stretch and a ‘growth stretch’. Following the notation of [37], we denote the scalar total and growthstretches by λ i and γ i ( i = s, φ, n ), respectively, so that λ s = α s γ s , (5a) λ φ = α φ γ φ , (5b) λ n = α n γ n , (5c)where the α i ( i = s, φ, n ) are the purely elastic stretches defined in equation (3).As defined in section 2.1, σ and ρ denote the arclength and radial distance in the unloaded state. That is, σ and ρ refer to the ‘virtual configuration’ that includes the part of the deformation due to growth, but notthe elastic response to the applied load. Recalling Σ as the arclength in the initial, ungrown and unloadedconfiguration, and denoting the associated radial distance by R , then the total and growth stretches in the e s and e φ directions may be expressed simply as λ s = ∂s∂Σ , (6a) γ s = ∂σ∂Σ , (6b) λ φ = rR , (6c) γ φ = ρR . (6d)We suppose that each material point in the sclera has some capacity for growth that could depend on a rangeof factors including age, position, and visual stimuli. We also assume that each region of the sclera can increasein size but cannot decrease, so that γ i ≥ i = s, φ, n ). It is known that the sclera becomes thinner duringaxial elongation, so, due to an absence of known mechanism, here we prohibit growth in scleral thickness andset γ n = 1, though alternative routes such as imposed mass conservation may easily be accommodated in thisframework. Thus, we have λ n = α n = h/H .In modeling the drivers of scleral growth, we are motivated by early experimental studies in which in-traocular pressure is associated with axial elongation and related myopia in chick embyros [40], rabbits [41],and humans [42], though the latter remains controversial. Hence, we assume that regions of the sclera growin response to elastic stress. Seeking a minimal phenomenological model, we assume that growth occurs torelieve the local strain, writing D γ s D t = η [ λ s − γ s ] + = ηγ s [ α s − + , (7a)D γ φ D t = η [ λ φ − γ φ ] + = ηγ φ [ α φ − + , (7b)where [ x ] + = max( x,
0) and η denotes a growth rate that is here assumed to be independent of direction,though may be generalised in future work. In equation (7), ‘D / D t ’ denotes a material derivative, holding Σ ,the arclength in the original configuration, fixed. These equations express that the relative rate of growth,˙ γ/γ , is driven by the elastic strains in the system. Since the intraocular pressure is held constant, with thesystem therefore always loaded, equilibrium is only achieved when η ( σ, t ) = 0, which we discuss later whenspecifying this growth rate.Since we anticipate that the final grown configuration will be far from the initial reference configurations, itis both preferable and computationally advantageous to express the dynamics in terms of the grown, unloadedconfiguration. Holding Σ constant, we write equation (7a) as ∂∂Σ (cid:18) D σ D t (cid:19) = η (cid:20) ∂s∂Σ − ∂σ∂Σ (cid:21) + . (8)Then, after multiplying through by ∂Σ/∂σ , applying the chain rule and integrating, we findD σ D t = Z σ η (cid:20) ∂s∂ ˆ σ − (cid:21) + dˆ σ , (9)noting that ∂Σ/∂σ ≥
0. As R is constant whilst Σ is constant, equation (7b) simply becomesD ρ D t = η [ r − ρ ] + , (10)as ∂R/∂ρ ≥
0. If η is constant, we recover the growth laws described in [39].5.4 Mechanical constitutive assumptionsThe timescale for scleral growth is on the order of years, so the equations of elastic equilibrium decouple fromthe growth laws and we consider equation (1), equation (2a), equation (3a) and equation (4) as a set of sevenordinary differential equations in the eleven unknowns r , z , s , θ , κ s , α s , q s , t s , t φ , m s and m φ , which eachdepend on σ . To close this sub-problem, we specify a constitutive relation for each of m s , m φ , t s and t φ ,thereafter viewing equations (1), (2a), (3a) and (4) as ordinary differential equations for r , z , θ , s , q s , α s and κ s . Firstly, we suppose that the curvature of the shell generates bending moments according to m s = m φ = E B ( κ s + κ φ ) , (11)where the bending modulus, E B , has units of force multiplied by length. This standard form of constitutivelaw is also used in [38, 39], though with a reference curvature. Here, we omit such a reference curvature, which,if constant, would have no impact on our model system, as we will see in section 2.6.Next, we view the sclera as an incompressible, hyperelastic, fibre-reinforced shell, accordingly assuming thatthe stress components t s and t φ are functions of the principal elastic stretches α s and α φ . A fibre-reinforcedformulation allows us to model the sclera’s ability to resist tension, as well as incorporating anisotropic effectsarising from fibre orientation, consistent with the substantial amount of collagen present in the sclera and theobservations of [15, 43].Fibre modeling in soft tissues is particularly challenging. When considered as a continuum, collagenoussoft-tissues show a certain amount of fibre distribution at each point [44–46]. In principle, this angular dis-tribution needs to be integrated at each point to obtain its overall mechanical contribution. However, whenthe distribution is sufficiently localized around different angles, the mechanical contribution of fibres can bemodeled using a finite number of reinforcing fibres [47]. The effect of a small fibre dispersion around theseparticular angles can also be taken into account, amounting to a renormalization of the elastic constants forboth the isotropic and anisotropic parameters [48, 49].Following our modeling philosophy of developing a tractable framework that captures the relevant me-chanical effects, we therefore model the anisotropic response of the sclera by introducing two families of fibres.However, since deformations have been assumed to be axisymmetric, the two fibre families must be equal instrength and opposite in alignment with respect to the main axes. Specifically, the two families of fibres makean angle ψ with e φ , where the angle ψ is a function of position. Denoting the fibre directions as a and b inthe initial reference configuration, we explicitly take a = sin ψ ( Σ ) e s + cos ψ ( Σ ) e φ , (12a) b = sin ψ ( Σ ) e s − cos ψ ( Σ ) e φ . (12b)This reduces to purely circumferential reinforcement when ψ = 0, in line with the observations of [46, 50],with dispersion about this configuration being modeled by small values of ψ . Inherent to this general setupis the issue of isotropy at the base of the sclera, where the surface intersects with its axis at Σ = σ = s = 0,which we resolve here by prescribing that ψ → π/ σ →
0, though this may also be addressed by omittinga small region of reinforcement close to this apex. Note that, here and throughout, fibre orientation is treatedas a material property, with ψ ( σ ( Σ, t ) , t ) = ψ ( Σ ) , (13)which we will later specify based on the observations of [15, 43]. We similarly treat H , the undeformed shellthickness, as a material property, in particular following [11], though we will also consider shells of initiallyuniform thickness for comparison.We now relate stress to strain via a strain-energy density, W , which, for simplicity, depends only on thefirst invariant of the isotropic strain, I , and the fibre stretch, I : I = α s + α φ + α n , (14a) I = α s sin ψ + α φ cos ψ . (14b)The principal stresses are then given by t s = 2 Hα n (cid:20)(cid:0) α s − α n (cid:1) ∂W∂I + 2 α s sin ψ ∂W∂I (cid:21) , (15a) t φ = 2 Hα n (cid:20)(cid:0) α φ − α n (cid:1) ∂W∂I + 2 α φ cos ψ ∂W∂I (cid:21) . (15b)6his is a consequence of standard shell theory assumptions concerning the shell’s thin aspect ratio [51],similar to the approach used, without fibre-reinforcement, in [38] and discussed with fibre-reinforcement in[52]. Further details are provided in appendix D.Several different strain-energy densities have been used to model sclera [15, 46, 53]. Keeping with our min-imal approach, here we adopt the simplest strain-energy density for an incompressible elastic fibre-reinforcedmaterial. Following many authors, we use a reinforced neo-Hookean strain-energy density of the form W = C ( I −
3) + D (cid:0) [ I − + (cid:1) , (16)where the first term represents the isotropic contribution to the stress and the second term describes the fibrereinforcement. Both C and D are material parameters with units of stress, and C can be identified with halfof the shear modulus in an isotropic neo-Hookean material. We take C within the range of values reportedin [54], whilst a range of values for D will be considered. In future work, both C and D could be allowed tovary with space and time, modeling for instance scleral hardening or softening with age or emmetropization,respectively. If the strains remain sufficiently low, this strain energy represents the dominant contribution ofmore sophisticated strain-energy functions that only depend on I and I , as it can be obtained as the first twoterms in a systematic expansion of a potential W = W ( I , I ). We also note that the presence of the invariant I is known to be important in shear problems, though is absent from this system due to axisymmetry.Substituting for W and noting that incompressibility ensures that α s α φ α n = 1, we find t s = 2 H (cid:18) C (cid:18) α s α φ − α s α φ ) (cid:19) + D [ I − + α s α φ sin ψ (cid:19) , (17a) t φ = 2 H (cid:18) C (cid:18) α φ α s − α s α φ ) (cid:19) + D [ I − + α φ α s cos ψ (cid:19) . (17b)2.5 Growth rateA plethora of stimuli and responses have been proposed that could combine in the human sclera to controlemmetropization. Our model allows us to easily simulate and compare these hypotheses, but in its firstexposition we restrict our attention to a simple scenario in order to illustrate the capabilities of the framework.Here, the local growth rate, η , represents the rate at which a region of the sclera will grow in order to relievethe local stretch. We decompose this growth rate into the product of an intrinsic capacity for growth, g c , anda stimulus response, g v . Explicitly, we write η ( σ, t ) = g c ( Σ ) g v ( σ, t ) , (18)where g c is a material property dependent only on Σ = Σ ( σ, t ). As we are not modeling the growth of thecornea, matching at the front of the eye requires g c → Σ → L , leading us to pose the phenomenologicalfunctional form g c ( Σ ) = η (cid:18) (cid:18) Γ − Σδ (cid:19)(cid:19) , (19)where the location and extent of the growing zone are governed by Γ and the decay away from this zone aswe move towards the anterior sclera is given by δ . The parameter η has units of inverse time and quantifiesthe maximum growth rate of the sclera. The stimulus dependence of the growth rate, g v , is defined in termsof the discrepancy between the current, deformed position of the sclera and some target surface, which willbe defined by the optical properties of the eye in its configuration at time t . Thus, g v is better understood inthe form g v ( σ, t ) = g v ( r ( σ, t ) , z ( σ, t ) , t ) , (20)recalling that s = s ( σ, t ).Adopting the hypothesised wavelength-dependent growth response of the sclera, we suppose that the retinacan detect the blurring of red and blue light that reaches it, and that this stimulates or inhibits a growthresponse. We model this simply, computing the best-focus surface for blue light as described in appendix A,then defining the growth response based on the location of the retina relative to this target surface. We thushave two cases to distinguish. If a material point on the retina is in front of the best-focus surface for bluelight, experiencing hyperopic defocus, then growth is triggered and its amplitude g v is set to be the distancebetween the material point and the closest point on the blue surface. Alternatively, if a point on the retina is7ehind the best-focus surface for blue light, experiencing myopic defocus, then growth is stopped and we set g v to zero. After simulation, the position of the retina is compared with the best-focus surface for red light todetermine if the eye is emmetropic. In particular, if the retina lies between the best-focus surfaces then we willterm the eye emmetropic, whilst we will describe it as myopic if the retina is behind the best-focus surface forred light and hyperopic if the retina lies in front of the best-focus surface for blue light.2.6 Model reduction and summaryFirst, we note that q s appears in the model only in combination with r , so we replace these terms with thenew variable Q = rq s . It is also convenient to write the equations of mechanical equilibrium in terms of thegrown, undeformed arclength σ , so that, making use of equation (3a), equations (1), (2a) and (4) become ∂r∂σ = α s cos θ , (21a) ∂z∂σ = − α s sin θ , (21b) ∂θ∂σ = α s κ s , (21c) ∂Q∂σ = α s [ r∆P − κ s rt s − κ φ rt φ ] , (21d) ∂t s ∂σ = α s r [ κ s Q + ( t φ − t s ) cos θ ] , , (21e) ∂m s ∂σ = α s r [( m φ − m s ) cos θ − Q ] . (21f)Substituting our constitutive assumptions for the bending moments into equation (21f), we find ∂κ s ∂σ = α s (cid:18) ( κ φ − κ s ) r cos θ − QrE b (cid:19) . (22)Equations (21a), (21c) to (21e) and (22) form a system of five ordinary differential equations in the fiveunknowns r , θ , α s , Q and κ s , each as a function of σ , noting that t s may be written as a function of α s .Following their solution, the variables z and s can be calculated using equation (3a) and equation (21b), i.e. z and s decouple. In the absence of experimental data to guide a choice of growth rate, η , we nondimensionalisetime by a typical timescale of growth, on the order of years, hereafter considering both time and growth rateto be dimensionless. Accordingly, we present simulations over a unit time interval.2.7 Initial and boundary conditionsThe ungrown unloaded reference configuration is parameterised by the arclength Σ ∈ [0 , L ], where Σ = L identifies the location where the sclera meets the cornea. The growing reference domain is thus σ ∈ [0 , σ ( L, t )].For simplicity, we assume that the original unloaded configuration has a spherical shape of radius R , so that R = R sin( Σ/R ). Thus, the initial conditions for equations (9) and (10) are simply given by ρ ( Σ,
0) = R sin (cid:18) ΣR (cid:19) , (23a) σ ( Σ,
0) =
Σ . (23b)We must also provide the initial scleral thickness and fibre orientation. In our simulations, we consider bothshells of uniform thickness and others with a more biologically realistic profile, with the sclera being thickestat the anterior pole, thinning towards the equator, then thickening towards the limbus. We give the explicitform of this non-uniform H in appendix B, in line with [11]. Based on observations in [15, 43], we consider afibre orientation that has greater reinforcement in the e s direction around the equatorial region of the sclera,but where fibres are predominantly aligned in the e φ direction at the limbus and near the posterior pole, againgiven explicitly in appendix B. 8t the back of the eye, where σ = 0, we impose r (0 , t ) = 0 , (24a) θ (0 , t ) = 0 , (24b) Q (0 , t ) = 0 , (24c)ensuring continuity of the shell and its slope, with the condition on Q being a consequence of requiring thenormal shear stress q s to be bounded as σ → σ = σ ( L, t ), we require the sclera to match smoothly to the cornea. The cornea is a deformable materialand therefore should stretch when the intraocular pressure ∆P is raised, suggesting that we should havesome pressure-dependent boundary condition. Instead of explicitly modeling the anterior eye, we model thecornea as a spherical material in its reference configuration. When pressure is increased, this surface deformsas detailed in appendix B, simply expanding as a spherical shell. The end points of the cornea are then usedas boundary conditions for the deforming sclera, which, as ∆P is held constant, are constant in time. Bymatching the surfaces, we have r ( σ ( L, t ) , t ) = r ⋆ , (25a) θ ( σ ( L, t ) , t ) = θ ⋆ , (25b)where the dependence of r ⋆ and θ ⋆ on the system parameters is explained in appendix B. After solving equa-tions (21a), (21c) to (21e) and (22) subject to the five boundary conditions in equation (24) and equation (25),we integrate equation (3a) and equation (21b) to obtain s and z . We set the centreline of the deformed cornealshell to be the origin of the frame, so that z ( σ ( L, t ) , t ) = z ⋆ , (26a) s (0 , t ) = 0 , (26b)where z ⋆ is defined in appendix B. Details of the implementation are given in appendix C, with typicalparameter values reported in table 1. η , withthe corresponding evolution of the retina relative to the best-focus surfaces for blue and red light shown infigure 3e. What is most notable is the qualitatively plausible shape attained by the sclera and retina, withthe local optically driven growth law apparently sufficient to produce realistic morphologies via the presentedminimal model. We also observe that the region of fastest growth migrates away from the posterior sclera, inline with the continued progress of the posterior retina towards the best-focus surface for blue light. Indeed,at later times the axial growth of the eye is being driven by growth in regions away from the posterior sclera,concentrated in a band nearing the scleral equator, with the rearmost portion of the retina having movedbehind the best-focus surface for blue light and therefore inhibiting local growth. In this case, this leads to amyopic eye, with the posterior retina at t = 1 having moved past the best-focus surface for red light.Figures 3f and 3g depict the elastic stretches in the e s and e φ directions, initially uniform and equal dueto the homogeneous spherical initial condition. The non-uniform growth of the sclera breaks this homogeneity,resulting in stretches that vary around the eye. Perhaps unsurprisingly, we note that regions of high α s areapproximately coincident with regions of reduced α φ , and vice versa.3.2 Effects of a non-uniform reference thicknessNow incorporating the non-uniform reference thickness of the sclera, as given explicitly in equation (27), weevaluate the differences in the morphology and elastic stretches between this and the previous uniform case offigure 3. The initial and final loaded configurations of the sclera are shown in figure 4a, where the solid blacklines denote the upper and lower surfaces of the sclera, whilst the dashed lines correspond to the case withuniform reference thickness. We observe only a small difference in scleral and retinal positioning due to the9 ) t = 0 . b ) t = 0 . c ) t = 0 . d ) t = 1 . Fig. 3
Numerical solutions of equations (21a), (21c) to (21e) and (22), showing growth of a uniform-thickness, fibre-freesclera. (a-d) The scleral shell throughout the growth and deformation, shown cut at select timepoints, shaded by the growthrate η . (e) The evolution of the retina over time, corresponding to the sclera shown in (a-d), with the target best-focussurface for blue light and the best-focus surface for red light also shown. Snapshots of the stretch in the (f) e s direction and(g) e φ direction at times t = 0 , . , . , . ,
1, with arrows indicating increasing time, noting the uniform initial stretchesdue to the initially uniform shell thickness. The regions of largest growth rate migrate away from the posterior eye overtime, with the posterior region rapidly nearing the target surface and thus experiencing slower growth. Despite this negativefeedback, the axial length of the eye increases past the best-focus surface for red light near the fovea, resulting here inmyopia. Parameter values used in this simulation are C = 100kPa, D = 0, H = 0 . E b = 18kPamm , η = 70mm − , Γ = 8mm and δ = 4mm.(a) (b) (c) Fig. 4
Growth of a fibre-free sclera with spatially varying unloaded thickness, with H as defined in equation (27), comparedto a sclera with uniform thickness. (a) The initial and final deformed sclera, where the solid and dashed lines are the surfacesof the sclera in the non-uniform and uniform cases, respectively, with the non-uniform case shaded. The stretches in the(b) e s and (c) e φ directions are shown as solid and dashed curves for the non-uniform and uniform sclera, respectively,with lighter curves corresponding to the initial condition and heavier curves corresponding to the final configuration. Thedifference in scleral reference structure appears to minimally impact the evolving shape of the eye, though differences inthe associated stretches are evident. In particular, the non-uniform reference configuration serves to reduce both of thesestretches near the posterior eye ( σ = 0). The parameters used in this simulation are ∆P = 2kPa, C = 100kPa, D = 0, E b = 18kPamm , η = 70mm − , Γ = 8mm and δ = 4mm. non-uniform reference thickness, though the retina is marginally shifted forwards compared to the uniformcase due to the increased scleral thickness in the posterior eye. However, the elastic stretches associated withthese final deformed configurations, shown in figure 4b and figure 4c, are more significantly altered, with themaximum stretches in the posterior sclera reduced in the variable reference thickness case. The variation in thestretches can be attributed directly to the non-uniform reference thickness, with thinner regions experiencingcomparatively larger stretches. 10 a) (b)(c) (d) Fig. 5
Effects of fibre-reinforcement on scleral growth. (a) The final deformed position of the sclera for various levels offibre reinforcement. Inset is a magnified view of the posterior sclera, which reveals that the strongest fibres result in anincreased axial length and reduced scleral angle at the back of the eye, though both changes are marginal. (b) The evolutionof the axial length of the eye, measured from z = 0 to the surface of the posterior retina, revealing slower and prolongedgrowth for the stronger fibres, though this effect is non-monotonic in fibre strength (cf D = 0 and D = 100). The stretchesin the (c) e s and (d) e φ directions typically decrease with increased fibre strength, though a complex relationship appearsaround σ = 21, where the stretch in the e φ direction increases with strong fibre reinforcement. The scleral thickness is givenby equation (27) and the fibre orientation is given by equation (28). The parameters used are ∆P = 2 kPa, C = 100 kPa, E b = 18 kPa mm , η = 70 mm − , Γ = 8 mm, δ = 4 mm, taking the non-uniform reference thickness H of equation (27).The fibre strength ranges from D =0 kPa to 400 kPa, as stated in the legend without units, though the D = 200 kPasimulation has been omitted from (a) and (b) for clarity. D , presenting a selection of grown deformedsclera in figure 5a. Immediately evident is the minimal effect that this fibre reinforcement has on the finalscleral morphology, with only small differences present. That being said, the highest degree of reinforcementconsidered here results in increased axial length and a reduced scleral angle at the posterior, as shown inset infigure 5a. The axial growth dynamics are illustrated in figure 5b, from which we again see the increased lengthof the most reinforced sclera, though resultant from a lower rate of growth over an extended period of timecompared to sclera with lower levels of reinforcement. Further, whilst the most reinforced shell has the greatestfinal length, the shell with only slight reinforcement exhibits a reduced length compared to the fibre-free shell.Thus, there is a non-monotonic response of the scleral morphology to the strength of fibre reinforcement. Theexistence of such a complex response of a scleral shell to fibre reinforcement is further illustrated in figure 5cand figure 5d, with the stretch in the e φ direction in the anterior sclera not following the overall trend of beingreduced by increased reinforcement.3.4 Impacts of growth lawsIn order to investigate the role of the intrinsic growth capacity, we consider a variety of parameter combinations( Γ, δ ), which parameterise the location of the growing region and the smoothness of the transition betweenregions of high and low growth capacity, respectively. Figure 6 shows significant changes in the morphologyof the grown deformed retina as we vary these parameters. At low values of δ , when the transition from lowto high growth capacity is rapid, we observe a marked change in retinal shape as the location of the growingzone is moved away from the posterior eye as Γ increases, with an upwards bump developing for Γ = 12 mm.As Γ is increased, we also note a change in the rate of growth, with the eye rapidly growing for high values11 = 4 Γ = 8 Γ = 12 δ = 2 δ = 6 Fig. 6
Effect of the growth zone on the dynamics and position of the retina, shown in the initial configuration, at t = 0 . t = 1, shaded by growth rate η . Varying both the location of the growing region, Γ , and the rate oftransition between low and high growth, δ , we observe a range of growth dynamics and final configurations of the retina. Asmoother transition between regions of fast and slow growth appears to result in a more spherical retina, whilst increasing Γ is associated with rapid initial growth rate, which is discernible from the location of the t = 0 . Γ do not always result inincreased axial length, as more clearly displayed in figure 7. Here, we have taken the non-uniform reference scleral thicknessof equation (27), considering fibre reinforcement as described in equation (28). The material parameters used here are ∆P = 2kPa, C = 100kPa, D = 100kPa and E b = 18kPamm , with η = 70 mm − . Both Γ and δ have units of millimetres. of Γ , though still attaining similar axial lengths at t = 1. This is shown explicitly in figure 7, in which theaxial growth rate can be seen to be strongly dependent on Γ when δ = 2 mm. Here, growth at the peripheralretina drives the axial progression at late stages of development, as seen earlier in figure 3. Here, we noteanother non-monotonic dependence of the ocular morphology on the parameters, that of final axial length onthe location of the growing region.Surprisingly, similar dependence on Γ is not present when considering δ = 6 mm, with the axial lengthof the eye at t = 1 approximately independent of Γ now that the intrinsic growth capacity transitions moregradually. Returning to figure 6, we also see that the final shape of the retina is largely unaffected by changesto Γ , despite the change in growth rate. Thus, the smoothness of the transition between regions of high and lowgrowth capacity appears dominant over the location of maximal growth and the associated rate of development.This suggests partial robustness of the retinal/scleral development process to the details of growth, furthersupported by the approximately consistent axial length seen throughout each of these simulations, with theexception of ( Γ, δ ) = (8 mm , δ = 6 mm, in starkcontrast to the varied shapes seen for δ = 2 mm. This suggests that the smoother transition of the growingregion in the former seemingly drives this more uniform growth and deformation, consistent with intuitionand distinct from the sharply varying growth dynamics present when δ = 2 mm.12 = 2 δ = 6 Fig. 7
Change in axial length of the eye over time for various regions of growth, corresponding to the configurations shownin figure 6. With the intrinsic capacity for growth specified as in equation (19), we see that increased Γ results in a fasterinitial rate of growth, though for high δ the maximum axial length is approximately independent of the location of maximalintrinsic growth, Γ . Indeed, the overall axial length at t = 1 appears to be robust to variations in the parameters, with allbut the ( Γ, δ ) = (12 ,
2) case being within a range of 3 mm. The material parameters used are as described in figure 6. Both Γ and δ have units of millimetres. This work has described and showcased an idealised model for the growth and development of the primarystructural component of the eye, the sclera. Under the assumption of morphoelasticity and axisymmetry, wehave derived a simple yet detailed model in which the growth mechanics of the thin scleral shell may be readilycoupled to external stimuli or material properties. Reducing to a simple system of five quasistatic ordinarydifferential equations, with the growth and elastic timescales separated by orders of magnitude, this flexibleframework may be solved with standard numerical methods, achievable without significant computational costor optimisation. This model is therefore well-suited to explorative studies of ocular development, sacrificingthe accuracy of geometrically refined models in favour of rapidly querying the fundamental principles that linkthe growth and deformation during emmetropization.In order to showcase the flexibility and utility of this approach, we have included and briefly explored theeffects of multiple mechanisms and features of the developing eye, focussing in particular on a hypothesiseddriver of growth. With the optical properties of the mid-growth eye either stimulating or inhibiting thegrowth of the model sclera based on the detection of hyperopic blur, we have seen that this hypothesiscan lead to qualitatively realistic ocular morphologies for typical and estimated parameter values. Havingsought throughout to impose the simplest plausible assumptions and constitutive laws for the geometry,mechanical properties and growth kinetics of the sclera, we have therefore seen that these are sufficient tophenomenologically capture the growth dynamics of the eye. In particular, the prescribed local growth lawwas sufficient to achieve an appropriate global response to external stimuli, forming eyes of a plausible sizeand shape for focused vision.Further, the observation of qualitative realism was generally found to be robust to changes in the detailsof the growth specification, though significant variation within this broad class was observed when modifyingthe intrinsic growth capacity of the sclera. This resulted in a range of configurations at maturity, most notablebeing the approximately spherical shapes observed when local growth varied more slowly over the sclera.Surprisingly, the inclusion of fibres had little effect in comparison to that of varying the growth specification.Indeed, the key property of the mature eye, the axial length, was largely unchanged by varying the degreeof reinforcement. However, the small observed changes exhibited a non-monotonic response to increases inreinforcement strength, suggesting the existence of a complex relationship between shell structure and growthdynamics that future study is expected to explore in detail.Our presented simulations posit a finite time interval in which the sclera can grow; thereafter, all growthstops. In practice, the growth rate may vary more smoothly with age and there may be stages in postnataldevelopmental when the eye experiences more rapid growth, for example. These non-uniform growth periodscould be considered via a simple extension of the framework developed in this work, and therefore could easilybe investigated. Additional readily realisable refinements could include coupling this model with the evolvingoptical properties of the front of the eye, with an interesting application being the evolution of best-focussurfaces during childhood and how the sclera grows in response to these changes. Another focus that meritsfurther theoretical investigation concerns the mechanism by which the sclera detects and responds to visual13 ig. 8
Computing the best-focus surfaces. A dense array of parallel, coherent rays (thin lines) is traced through the anterioroptics (cornea, pupil, lens – medium lines) and the phase of the wavefront emerging through the posterior surface of thelens is fitted to Zernike polynomial functions. Bottom – Left: Example of a lens-emerging waveform. A Kirchhoff integralis applied to further propagate the wavefront and compute its modulation transfer function (MTF) on small surfacesperpendicular to the central ray. Bottom – Right: Example MTFs along probed surfaces. Surface C maximizes the Strehlratio and hence corresponds to the best-focus surface for this wavelength. Computing the best-focus distance for a rangeof incidence angles (0-40 degrees), assuming axial symmetry, we reconstruct the best-focus surfaces for two wavelengths:400 nm (blue) and 600 nm (red). information, for example the effects of individual variation in photoreceptor topography and its link to theaxial length of the eye [55]. If chromatic effects are demonstrated to be significant in emmetropization, onecan envisage the design of corrective eyewear, guided by mathematical modeling, that specifically tailors theshape of image surfaces to moderate the growth of the eye.In summary, we have presented a simple morphoelastic model of the complex multifaceted process ofemmetropization. The model provides a step towards improving understanding of this developmental process,demonstrating that local growth laws can lead to qualitatively realistic morphologies of the elastically deformedeye, whilst enabling simple yet detailed future explorations of varied hypotheses for ocular development.
Acknowledgements
The research leading to these results has received funding from the European Union Seventh Frame-work Programme (FP7/2007-2013) under grant agreement no. 309962 (HydroZONES). BJW is supported by the UKEngineering and Physical Sciences Research Council (EPSRC), Grant No. EP/N509711/1.
The computer code used and generated in this work is freely available from https://gitlab.com/bjwalker/morphoelastic-eye.git
A Optical calculations
The optical components of the anterior eye, the cornea, anterior chamber, lens and vitreous chamber, focus the lightwavefronts that are incident on the eye on a fictitious curved surface near the retina, which we term the best-focus surface.The position and shape of this surface are dependent on the wavelength of the incident light due to chromatic aberrationsin the light focusing components. Modeling the geometrical and optical properties of the anterior eye as in [56], a raytracingalgorithm was employed in order to compute the individual surfaces of best focus for red and blue incident light wavefronts,exemplified in figure 8. Dense arrays of parallel coherent rays were traced through the anterior optics, with the phase ofthe wavefront emerging on the posterior surface of the lens fitted to Zernike polynomial functions. These are propagatedvia a Kirchhoff integral and the Strehl ratio is computed on various test surfaces perpendicular to the central ray. Thesurface corresponding to the maximum Strehl ratio represents the best-focus surface, which is constructed for incident anglesbetween 0 and 40 degrees, appealing to assumed axisymmetry. This approach may be readily extended to include the effectsof additional or non-uniform lenses, enabling the modeling of corrective lenses and their effects on ocular development, forexample. a) (b) Fig. 9
Reference non-uniform scleral thickness and fibre orientation, shown in (a) and (b), respectively, from the works of[11] and [15].
B Initial and boundary conditions
When considering a non-uniform scleral thickness, following [11] we prescribe H ( Σ ) = (cid:26) . − . (cid:0) Σ − . πR . (cid:1) , x ∈ [0 , . πR ] , .
475 + 0 .
025 tanh( Σ − . πR ) , x ∈ (0 . πR , L ] , (27)as shown in figure 9 alongside the initial fibre orientation, prescribed as ψ ( Σ ) π = 0 .
25 + , x ∈ [0 , . πR ] , .
06 cos (cid:0) Σ − . πR . R (cid:1) − . , x ∈ (0 . πR , . πR ] , . − .
17 cos (cid:0) Σ − . πR . R (cid:1) , x ∈ (0 . πR , . πR ] , .
22 cos (cid:0) π ( Σ − . πR ) L − . πR (cid:1) , x ∈ (0 . πR , L ] , (28)following the observations of [15, 43].At the anterior point of the sclera, we match the scleral displacement to the inflation of a thin, spherically symmetric,non-growing shell of uniform thickness with no fibres, minimally modeling the cornea. Firstly, for this simple shell, we seethat α cs = α cφ , so for notational convenience we denote the stretch simply by α , where the superscript on all other variablesdenotes that we are considering the cornea. Since the corneal reference configuration is spherical, we have R c = R c sin (cid:18) Σ c R c (cid:19) , (29)for Σ c ∈ [0 , πR c ], where R c is the radius of the sphere. The position of a point on the inflated sphere is thus given by r c = αR c sin (cid:18) Σ c R c (cid:19) , (30a) z c = αR c cos (cid:18) Σ c R c (cid:19) + B , (30b)where B is a constant of integration. The constraint of spherical symmetry ensures there is no normal shear force, Q c = 0,so that the shell deforms as a membrane. The solution to this problem is presented in [57], where it is shown that ∆P = 4 C c H c αR c (cid:16) − α (cid:17) , (31)where C c is the neo-Hookean constant, H c is the undeformed thickness and R c is the undeformed radius of the shell. Wetake C c , H c and R c to have values based on the mechanics of the cornea and we calculate α for the required pressuredifference numerically, restricting α ∈ (1 , (1 / ) due to the non-injective relation between α and ∆P . In particular, theupper limit here is the α value corresponding to the maximum of 1 /α − /α for α >
1, it placing a bound on the maximumpressure difference that we can consider, though we don’t vary the intraocular pressure in this work. Finally, evaluatingthe shape of the cornea at the point of attached to be sclera, we find the boundary conditions for the scleral shell to be r ⋆ = αR c sin (cid:18) LR c (cid:19) , (32a) θ ⋆ = LR c , (32b) z ⋆ = αR c (cid:18) (cid:18) LR c (cid:19)(cid:19) . (32c) ame Value Source R
10 mm [58] L . πR mm [58] H .
65 mm Based on equation (27) ∆P C
100 kPa [54] D
100 kPa to 400 kPa Estimated E b
18 kPa mm [59] η
70 mm − Estimated δ Γ Table 1
Typical parameter values used in simulations, unless otherwise specified. R and L represent the size of a typicalsclera towards the end of the rapid growth phase of development, at the beginning of the time interval that we model. Theestimate for C is based on values measured in monkey sclera, albeit with a slightly different constitutive law. The bendingstiffness is estimated from C and H under the assumption of incompressibility. The corneal parameters R c , H c and C c arechosen to match the scleral properties. C Implementation
The governing equations presented in section 2.6 have a singularity when r = 0, so we solve the system numerically onthe truncated domain σ ∈ [ σ ( ε, t ) , σ ( L, t )] for 0 < ε ≪
1. By expanding the variables r , θ , α s , κ s and Q around σ = 0and evaluating at σ = σ ( ε, t ), following [39], then substituting the expansions into (21 a,c,d,e ) and equation (22) subjectto equation (24), it becomes clear that the singularity is removable for compatible initial fibre directions. Indeed, isotropyis required as σ → ψ → π/ σ → σ → ψ → π/ σ →
0, subject to which the stress is finite and theboundary conditions on r and θ can be replaced by the notationally cumbersome r ( ε ) = σ ( ε ) α s ( σ ( ε )) , (33a) θ ( ε ) = σ ( ε ) α s ( σ ( ε )) κ s ( σ ( ε )) , (33b)where we have suppressed the t -dependence of all quantities here for brevity. Now considering Q , we further manipulateequation (21) to admit the first integral Q cos θ + rt s sin θ − r ∆P A, (34)where A is a constant. Since we require solutions that pass through r = 0 with θ = π/
2, we find A = 0. Thus, evaluatingequation (34) at σ = σ ( ε, t ) provides the analogous truncated boundary condition for Q . Note that whilst it is possible touse equation (34) to eliminate Q from equations (21a), (21c) to (21e) and (22), preliminary numerical simulations suggestedthat it is easier to solve the five ordinary differential equations than the reduced system. Hence, we retain Q in the governingequations and use equation (34) as a check on the numerical solutions.We utilise MATLAB’s inbuilt adaptive boundary problem solver bvp4c to solve equations equations (21a), (21c) to (21e)and (22) subject to the boundary conditions truncated boundary conditions. The initial conditions are provided on a regulargrid for Σ ∈ [0 , L ] and growth is approximated with an explicit Euler scheme for equations (9) and (10). For each simulation,we ensure that the solution has converged with respect to our choices of grid size, timestep, truncation point and errortolerances in the solver. For example, the simulations in figure 3 were rerun on a refined spatial grid, with a smallertimestep, with a lower error tolerance in the bvp4c solver, and with a reduced truncation value ε . The largest relative errorsin the variables κ s , α s , r and θ at t = 1 in this refined simulation are 1 . × − , 2 . × − , 9 . × − and 7 . × − respectively, well below practical tolerance. Typical parameter values for the simulations in this work are given in table 1. D Shell stresses
We briefly discuss the standard rationale for the functional form of the stresses specified by equations (17a) and (17b). Ina fully three-dimensional elastic body, the Cauchy stress, σ , is σ = − p I + 2 F ∂W∂ C F T , (35)where F is the deformation gradient, C = F T F is the right Cauchy-Green tensor, p is the hydrostatic contribution to thestress associated with enforcing incompressibility, and W is the strain-energy function. For a detailed account, we directthe interested reader to, for example, the work of [52]. If we suppose that W = W ( I , I , I ) where I = tr C , (36a) I = a T Ca , (36b) I = b T Cb , (36c) o that I and I represent the stretch of fibres that lie in the directions a and b , then σ = − p I + 2 F (cid:16) ∂W∂I I + ∂W∂I a ⊗ a + ∂W∂I b ⊗ b (cid:17) F T . (37)Furthermore, if we select the basis used for our scleral model so that F is diagonal with entries α s , α φ , α n , define the fibredirections as in equations (12a), (12b) and (13), so that I = I , and finally require W ( I , I , I ) = W ( I , I , I ), then theoff-diagonal terms in a ⊗ a and b ⊗ b cancel. Thus, the only non-zero components in equation (37) are σ ss = − p + 2 α s ∂W∂I + 4 α s sin ψ ∂W∂I , (38a) σ φφ = − p + 2 α φ ∂W∂I + 4 α φ cos ψ ∂W∂I , (38b) σ nn = − p + 2 α n ∂W∂I . (38c)The shell’s thin geometry can be exploited as discussed in the context of membranes in [60]. We apply the key results inour shell model by working with resultant stresses of the form t s = α n Hσ ss , (39a) t φ = α n Hσ φφ (39b)and setting σ nn = 0, often termed the ‘membrane assumption’. This is akin to noting that the curved shell is so thin thatload across the surface due to the intraocular pressure is supported by in-shell tension, as opposed to stress across theshell thickness. The membrane assumption enables the elimination of the hydrostatic pressure, p , and our incompressibilityassumption, α n = 1 /α s α φ , gives principal in-shell stress resultants of the form equation (14). References
1. Morgan, I.G., Ohno-Matsui, K., Saw, S.M.: Myopia. Lancet (9827), 1739–1748 (2012)2. Spillmann, L.: Stopping the rise of myopia in Asia. Graefe’s Archive for Clinical and Ex-perimental Ophthalmology (5), 943–959 (2020). DOI 10.1007/s00417-019-04555-0. URL http://link.springer.com/10.1007/s00417-019-04555-0
3. Sherwin, J.C., Mackey, D.A.: Update on the epidemiology and genetics of myopic refractive error. Expert Review ofOphthalmology (1), 63–87 (2013)4. Dolgin, E.: The myopia boom. Nature (7543), 276–278 (2015)5. Holden, B.A., Jong, M., Davis, S., Wilson, D., Fricke, T., Resnikoff, S.: Nearly 1 billion myopes at risk of myopia-relatedsight-threatening conditions by 2050 - time to act now. Clinical and Experimental Optometry (6), 491–493 (2015).DOI 10.1111/cxo.12339. URL http://doi.wiley.com/10.1111/cxo.12339
6. Holden, B.A., Fricke, T.R., Wilson, D.A., Jong, M., Naidoo, K.S., Sankaridurg, P., Wong, T.Y., Naduvi-lath, T.J., Resnikoff, S.: Global Prevalence of Myopia and High Myopia and Temporal Trends from2000 through 2050. Ophthalmology (5), 1036–1042 (2016). DOI 10.1016/j.ophtha.2016.01.006. URL https://linkinghub.elsevier.com/retrieve/pii/S0161642016000257
7. Morgan, R., Speakman, J., Grimshaw, S.: Inuit myopia: an environmentally induced “epidemic”? Can. Med. Assoc. J. (5), 575 (1975)8. Manjunath, V., Enyedi, L.: Pediatric myopic progression treatments: Science, sham, and promise. Current Ophthal-mology Reports (4), 150–157 (2014)9. Russo, A., Semeraro, F., Romano, M.R., Mastropasqua, R., Dell’Omo, R., Costagliola, C.: Myopia onset and progres-sion: can it be prevented? Int. Ophthalmol. (3), 693–705 (2014)10. Wu, P.C., Chuang, M.N., Choi, J., Chen, H., Wu, G., Ohno-Matsui, K., Jonas, J.B., Cheung, C.M.G.: Update in myopiaand treatment strategy of atropine use in myopia control. Eye (1), 3–13 (2019). DOI 10.1038/s41433-018-0139-711. Fatt, I., Weissman, B.: Physiology of the eye: an introduction to the vegetative functions, 2nd edn. ButterworthsLondon (1992)12. Karimi, A., Razaghi, R., Navidbakhsh, M., Sera, T., Kudo, S.: Mechanical Properties of the Human Sclera Under Vari-ous Strain Rates: Elastic, Hyperelastic, and Viscoelastic Models. Journal of Biomaterials and Tissue Engineering (8),686–695 (2017). DOI 10.1166/jbt.2017.1609. URL
13. Romano, M.R., Romano, V., Pandolfi, A., Costagliola, C., Angelillo, M.: On the use of uniaxial tests on the sclerato understand the difference between emmetropic and highly myopic eyes. Meccanica (3), 603–612 (2017). DOI10.1007/s11012-016-0416-0. URL http://link.springer.com/10.1007/s11012-016-0416-0
14. Grytz, R., Girkin, C.A., Libertiaux, V., Downs, J.C.: Perspectives on biomechanical growth and remodeling mechanismsin glaucoma. Mech. Res. Commun. , 92–106 (2012)15. Grytz, R., Fazio, M.A., Girard, M.J., Libertiaux, V., Bruno, L., Gardiner, S., Girkin, C.A., Downs, J.C.: Materialproperties of the posterior human sclera. J. Mech. Behav. Biomed. , 602–617 (2013)16. Grytz, R., El Hamdaoui, M.: Multi-scale Modeling of Vision-Guided Remodeling and Age-Dependent Growth of the Tree Shrew Sclera During Eye Development and Lens-Induced My-opia. Journal of Elasticity (1-2), 171–195 (2017). DOI 10.1007/s10659-016-9603-4. URL http://dx.doi.org/10.1007/s10659-016-9603-4http://link.springer.com/10.1007/s10659-016-9603-4
17. Simonini, I., Pandolfi, A.: Customized Finite Element Modelling of the Human Cornea. PLOS ONE (6), e0130426(2015). DOI 10.1371/journal.pone.0130426. URL https://dx.plos.org/10.1371/journal.pone.0130426
8. S´anchez, P., Moutsouris, K., Pandolfi, A.: Biomechanical and optical behavior of human corneas before and afterphotorefractive keratectomy. Journal of Cataract & Refractive Surgery (6), 905–917 (2014)19. Pandolfi, A.: Cornea modelling. Eye and Vision (1), 1–15 (2020). DOI 10.1186/s40662-019-0166-x20. Humphrey, J.D.: Cardiovascular solid mechanics. Cells, tissues, and organs. Springer Verlag, New York (2002)21. Cowin, S.C.: Tissue growth and remodeling. Annu. Rev. Biomed. Eng. , 77–107 (2004)22. Gasser, T.C., Ogden, R.W., Holzapfel, G.A.: Hyperelastic modelling of arterial layers with distributed collagen fibreorientations. J. R. Soc. Interface (6), 15–35 (2006)23. Kuhl, E., Maas, R., Himpel, G., Menzel, A.: Computational modeling of arterial wall growth. Biomech. Model. Mechan. (5), 321–331 (2007)24. Taber, L.A.: Biomechanics of growth, remodeling and morphogenesis. Appl. Mech. Rev. , 487–545 (1995)25. Moulton, D.E., Goriely, A.: Possible role of differential growth in airway wall remodeling in asthma. J. Appl. Physiol. (4), 1003–1012 (2011)26. Budday, S., Steinmann, P., Goriely, A., Kuhl, E.: Size and curvature regulate pattern selection in the mammalianbrain. Extreme Mech. Lett. , 193–198 (2015)27. Diether, S., Schaeffel, F.: Local changes in eye growth induced by imposed local refractive error despite active accom-modation. Vision Res. (6), 659–668 (1997)28. Wallman, J., Gottlieb, M.D., Rajaram, V., Fugate-Wentzek, L.A.: Local retinal regions control local eye growth andmyopia. Science (4810), 73–77 (1987)29. Foulds, W.S., Barathi, V.A., Luu, C.D.: Progressive myopia or hyperopia can be induced in chicks and reversed bymanipulation of the chromaticity of ambient light. Invest. Ophth. Vis. Sci. (13), 8004–8012 (2013)30. Kr¨oger, R., Wagner, H.J.: The eye of the blue acara (aequidens pulcher, cichlidae) grows to compensate for defocusdue to chromatic aberration. J. Comp. Physiol. A (6), 837–842 (1996)31. Liu, R., Qian, Y.F., He, J.C., Hu, M., Zhou, X.T., Dai, J.H., Qu, X.M., Chu, R.Y.: Effects of different monochromaticlights on refractive development and eye growth in guinea pigs. Exp. Eye Res. (6), 447–453 (2011)32. Werblin, F., Roska, B.: The movies in our eyes. Sci. Am. (4), 72–79 (2007)33. Oyster, C.W.: The human eye. Sinauer Associates (1999)34. Rodriguez, E.K., Hoger, A., McCulloch, A.D.: Stress-dependent finite growth in soft elastic tissues. J. Biomech. (4),455–467 (1994)35. Ambrosi, D., Ateshian, G., Arruda, E., Cowin, S., Dumais, J., Goriely, A., Holzapfel, G.A., Humphrey, J., Kemkemer,R., Kuhl, E., et al.: Perspectives on biological growth and remodeling. J. Mech. Phys. Solids (4), 863–883 (2011)36. Kuhl, E.: Growing matter: A review of growth in living systems. J. Mech. Behav. Biomed. , 529–543 (2014)37. Goriely, A.: The Mathematics and Mechanics of Biological Growth. Interdisciplinary Applied Mathematics. SpringerNew York (2017). URL https://books.google.co.uk/books?id=rgImDwAAQBAJ
38. Tongen, A., Goriely, A., Tabor, M.: Biomechanical model for appressorial design in magnaporthe grisea. J Theor. Biol. (1), 1–8. [We note there is a sign error in equation (8) corrected in Woolley et al. (2014)] (2006)39. Woolley, T.E., Gaffney, E.A., Oliver, J.M., Baker, R.E., Waters, S.L., Goriely, A.: Cellular blebs: pressure-driven,axisymmetric, membrane protrusions. Biomech. Model. Mechan. (2), 463–476 (2014)40. Coulombre, A.J.: The role of intraocular pressure in the development of the chick eye. i. control of eye size. Journal ofExperimental Zoology (2), 211–225 (1956)41. Maurice, D., Mushin, A.: Production of myopia in rabbits by raised body-temperature and increased intraocularpressure. The Lancet (7474), 1160–1162 (1966)42. Quinn, G.E., Berlin, J.A., Young, T.L., Ziylan, S., Stone, R.A.: Association of intraocular pressure and myopia inchildren. Ophthalmology (2), 180–185 (1995)43. Girard, M.J., Dahlmann-Noor, A., Rayapureddi, S., Bechara, J.A., Bertin, B.M., Jones, H., Albon, J., Khaw, P.T.,Ethier, C.R.: Quantitative mapping of scleral fiber orientation in normal rat eyes. Investigative Ophthalmology &Visual Science (13), 9684–9693 (2011)44. Jones, H., Girard, M., White, N., Fautsch, M.P., Morgan, J., Ethier, C., Albon, J.: Quantitative analysis of three-dimensional fibrillar collagen microstructure within the normal, aged and glaucomatous human optic nerve head.Journal of The Royal Society Interface (106), 20150066 (2015)45. Gouget, C.L., Girard, M.J., Ethier, C.R.: A constrained von mises distribution to describe fiber organization in thinsoft tissues. Biomechanics and modeling in mechanobiology (3-4), 475–482 (2012)46. Coudrillier, B., Pijanka, J.K., Jefferys, J.L., Goel, A., Quigley, H.A., Boote, C., Nguyen, T.D.: Glaucoma-relatedchanges in the mechanical properties and collagen micro-architecture of the human sclera. PloS one (7), e0131396(2015)47. Spencer, A.J.M.: Deformations of Fibre-Reinforced Materials. Oxford (1972)48. Holzapfel, G.A., Ogden, R.W.: Constitutive modelling of arteries. Proc. R. Soc. A (2118), 1551–1597 (2010)49. Melnik, A.V., Da Rocha, H.B., Goriely, A.: On the modeling of fiber dispersion in fiber-reinforced elastic materials.Int. J. Non Linear Mech. (2015)50. Zhang, L., Albon, J., Jones, H., Gouget, C.L., Ethier, C.R., Goh, J.C., Girard, M.J.: Collagen microstructural factorsinfluencing optic nerve head biomechanicscollagen microstructural factors. Investigative ophthalmology & visual science (3), 2031–2042 (2015)51. Green, A.E., Adkins, J.E.: Large elastic deformations. Oxford University Press (1970)52. Holzapfel, G.A., Ogden, R.W.: Constitutive modelling of arteries. P. R. Soc. A. (2118), 1551–1597 (2010)53. Coudrillier, B., Boote, C., Quigley, H.A., Nguyen, T.D.: Scleral anisotropy and its effects on the mechanical responseof the optic nerve head. Biomech Model Mechan (5), 941–963 (2013)54. Girard, M., Downs, J.C., Bottlang, M., Burgoyne, C.F., Suh, J.F.: Peripapillary and posterior scleral mechanics, partII-experimental and inverse finite element characterization. J. Biomech. Eng. (5), 051012 (2009)55. Wang, Y., Bensaid, N., Tiruveedhula, P., Ma, J., Ravikumar, S., Roorda, A.: Human foveal cone photoreceptor topog-raphy and its dependence on eye length. eLife , e47148 (2019). URL https://doi.org/10.7554/eLife.47148
56. Atchison, D.A.: Optical models for human myopic eyes. Vision Res. (14), 2236–2250 (2006)57. Adkins, J.E., Rivlin, R.S.: Large elastic deformations of isotropic materials. IX. The deformation of thin shells. Philos.T. R. Soc. S.-A. (888), 505–531 (1952)
8. Gordon, R.A., Donzis, P.B.: Refractive development of the human eye. AMA Arch. Ophthalmol. (6), 785–789(1985)59. Howell, P., Kozyreff, G., Ockendon, J.: Applied solid mechanics. Cambridge University Press (2009)60. Haughton, D.: Elastic membranes. London Mathematical Society Lecture Note Series pp. 233–267 (2001)(6), 785–789(1985)59. Howell, P., Kozyreff, G., Ockendon, J.: Applied solid mechanics. Cambridge University Press (2009)60. Haughton, D.: Elastic membranes. London Mathematical Society Lecture Note Series pp. 233–267 (2001)