The importance of geometry in the corneal micropocket assay

The corneal micropocket assay is an experimental protocol for studying vessel network formation, or neovascularization, in vivo. The assay is attractive due to the ease with which the developing vessel network can be observed in the same animal over time. However, it is difficultto relate observations of neovascularization in the hemispherical geometry of the cornea to predictions from mathematical (in silico) models and observations from in vitro models, which mostly adopt planar geometries. Here, we develop a three-dimensional, off-lattice mathematicalmodel of neovascularization in the cornea, with a new repulsion model for the interaction of vessels with tissue boundaries. The model allows us to: i) investigate how the hemispherical geometry of the cornea influences vessel network formation and ii) identify the metrics of neovascularization that are most robust to geometrical differences between typical in silico, in vitro and in vivo tissue domains. We predict that two-dimensional geometries lead to increased vessel ‘merging’, or anastomosis, relative to three-dimensional versions, with differences being most pronounced when circular geometries are adopted. Cuboidal geometries lead to patterns of neovascularization that most closely mimic the in vivo case. The ‘distance of the vascular front to the limbus’ is found to be a robust metric of neovascularization across the studied tissue domains, while densities of tip cells and vessels and ‘vascularized fraction’ are more sensitive to domain choice. Given the widespread adoption and attractive simplicity of planar tissue domains, the differences identified in the present study should prove useful in relating the results of previous and future in silico and in vitro studies of neovascularization to in vivo observations in the cornea. Author summary Neovascularization, or the formation of new blood vessels, is an important process in development, wound healing and cancer. The corneal micropocket assay is used to better understand the process and, in the case of cancer, how it can be controlled with drug therapies for improved patient outcomes. In the assay, the hemispherical shape of the cornea can influence the way the vessel network forms. This makes it difficult to directly compare results from experiments with the predictions of mathematical models or cell culture experiments, which are typically performed on flat substrates or planar matrices. In this study, we use mathematical modeling to investigate how the hemispherical shape of the cornea affects vessel formation and to identify how sensitive different measurements of neovascularization are to geometry.

typically performed on flat substrates or planar matrices. In this study, we use mathematical modeling to investigate how the hemispherical shape of the cornea affects vessel formation and to identify how sensitive different measurements of neovascularization are to geometry. Introduction 1 Neovascularization, or new blood vessel formation, is an important process in development, 2 wound healing, cancer and other diseases. The corneal micropocket assay, shown in Fig 1, is 3 widely used for studying neovascularization in vivo [1][2][3]. The assay involves the implantation 4 of a pellet containing pro-angiogenic compounds into the cornea of a small rodent and 5 observation of the resulting vessel formation over time. An attractive feature of the assay is the 6 relative ease with which neovascularization can be observed in the same animal at multiple 7 time-points. In rodents, the cornea is hemispherical, with a thickness on the order of 10 to 20 8 vessel diameters [4]. The pellet is not usually located on the pole of the cornea; it is shifted 9 slightly toward the base or limbus. It is evident from experimental images that the geometrical 10 configuration of the cornea-pellet system influences neovascularization patterns [1]. As shown 11 schematically in Fig 1C, new vessel formation is often focused in the region where the distance 12 between the pellet and limbus is smallest, with vessels tending to grow toward the pellet. 13 Fig 1. Schematics of the corneal micropocket assay and typical in silico and in vitro geometries. A) An image of the micropocket assay from Connor et al. [5]. B) The placement of a pellet of radius r p into the cornea at height h from the limbus, where h is the distance along the shown y or z axis. C) The resulting neovascularization, quantified by distance d of the front to the limbus. D) A typical planar model of the cornea-pellet system, with the pellet width w set equal to cornea width W for one-dimensional mathematical models. E) A typical circular model.
The attractive features of the micropocket assay make it a good candidate for comparing 14 vessel network formation with mathematical (in silico) and in vitro models. However, it is 15 important to understand how the geometry of the cornea-pellet system affects vessel network 16 formation relative to the planar tissue domains used in typical two-and three-dimensional (2D, 17 3D) cultures and mathematical models, a selection of which are shown in Figs 1D and E. 18 As reviewed in Jackson and Zheng [6], many mathematical models of neovascularization have 19 been motivated by the micropocket assay, providing valuable insights into the process [5,[7][8][9][10][11][12][13][14]. 20 To date, these models have adopted either one-dimensional (1D) or 2D representations of the 21 cornea-pellet system, with the former allowing efficient, continuum modeling of the developing 22 vessel network by the solution of systems of partial differential equations (PDEs) [5,10,11]. 23 However, the in vivo configuration of a cylindrical pellet placed 'off-center' in a hemispherical 24 volume is most naturally described using a 3D representation. While there are many 3D 25 mathematical models of sprouting angiogenesis [15][16][17][18][19], none have focused on the particular 26 geometry of the cornea-pellet system. Doing so brings the additional challenges of needing to 27 use a relatively large simulation domain and accounting for the interaction of vessels with 28 curved tissue boundaries at the epithelial and endothelial surfaces of the cornea. 29 In the present study, we develop a discrete, 3D, off-lattice mathematical model of 30 neovascularization in the cornea-pellet system, focusing on emulating the in vivo configuration. 31 We use a simplified treatment of the underlying biology, focusing instead on how the adoption 32 of different in vitro and in silico geometries, including planar 2D and 3D cultures in rectangular 33 or circular configurations, affects vessel network formation relative to the in vivo case. The  [20], recently developed by the 38 authors. 39 As part of this comparison we also identify metrics of neovascularization with high and low 40 sensitivity to the choice of tissue domain. As a result, we can predict which in vitro and in Neovascularization is simulated in seven different tissue domains, representing those typical of 48 in silico and in vitro modeling studies and the in vivo assay. A phenomenological model of 49 sprouting angiogenesis is adopted, motivated by several previous studies [6,7,14,21], but 50 extended to 3D. A 'soft-contact' model is also introduced to account for interactions of 51 migrating vessels with the cornea boundaries. It is assumed that the pellet contains a single 52 pro-angiogenic compound, Vascular Endothelial Growth Factor 165 (denoted VEGF in the 53 present study). Two situations are considered, one where the concentration of VEGF in the 54 tissue domain is described by a time-independent, spatially varying field where VEGF levels 55 decrease linearly from the pellet, and the other where VEGF dynamics are explicitly modeled, 56 following approaches in Tong and Yuan [7] and Connor et al. [5].

57
The simulated tissue domains, angiogenesis model and VEGF dynamics model are described 58 in this section, along with details of the numerical schemes used for their solution. used when presenting results. The pellet radius r p = 200 µm is based on data from Connor et 62 al. [5], but reduced from their value of 300 µm to facilitate placement in the cornea. In 3D 63 simulations, pellets are assumed to have a thickness of T p = 40 µm and are situated mid-way 64 between the epithelial and endothelial sides of the cornea. In 2D simulations, the cornea 65 thickness is neglected, while in 3D a value of T = 100 µm is used [4]. The cornea radius is fixed 66 at R = 1300 µm, which is a suitable value for mouse [4]. The 'Hemisphere' geometry is formed 67 by a 360°revolution of a circular arc of radius R and angle 90°about the polar axis, followed by 68 an extrusion through a distance T along the inward normal to the revolved surface, giving a 3D 69 volume. The cylindrical pellet is placed inside this volume and is completely enclosed by it. For 70 the Hemisphere, the pellet height h above the limbus is the distance as projected into 2D, as 71 would be typically measured in experimental images, rather than the distance along the 72 geodesic from the limbus to the pellet. 73 Fig 2. 3D renderings of the studied tissue domains, including dimensions. 3D renderings of the studied simulation domains, including adopted abbreviations. Domain dimensions are included in µm, following the notation in Fig 1, and T is the cornea thickness.
All simulations begin with a single blood vessel positioned a small height ε = 100 µm above 74 the base of the cornea and mid-way between the epithelial and endothelial sides of the cornea. 75 The vessel occupies the entire width (or circumference) of the domain at that position. Vessels 76 are discretized into straight line segments joined by 'nodes', shown schematically in Fig 3. Each 77 segment represents an arbitrary number of endothelial cells. Initially, nodes are placed 78 1 E L = 20.0 µm apart, where E L represents the number of endothelial cells per vessel length [22]. 79

80
A simple, phenomenological model of sprouting angiogenesis is used. The average rate of sprout 81 formation at a node located at x is [21]: where P max is the rate of sprouting per cell, c(x, t) is the VEGF concentration at location x 83 and time t, c 50 is the VEGF concentration at which the rate of sprouting is half-maximal and 84 L s is the averaged length of the two line segments joined to the node. In addition, a simple 85 description of lateral inhibition is used, with P = 0 within a distance 1 E L of a node that has 86 already been selected for sprouting. Simulations are discretized in time using a fixed step of 87 ∆t = 1.0 h. In each time step a random number z ∈ [0, 1] is chosen from a uniform distribution 88 at each node and a sprout forms if z < P ∆t. A different random number is generated at each 89 node.

90
After sprouting, migrating tips, illustrated in Fig 3, are assumed to move at constant speed 91 s = 10 µm h −1 . This speed is chosen so that the average time for the vascular front to reach the 92 pellet is on the order of 4 days, which is consistent with experimental observations [1]. A 93 persistent, off-lattice random walk is used to describe the migration of tips through the 94 extracellular matrix of the stroma [6,7,14]. The migration direction m is given by: where χ is a dimensionless weighting parameter controlling chemotactic sensitivity, m p is a unit 96 vector in the persistence direction and ∇c is the gradient of the VEGF concentration. The 97 random persistence direction m p is obtained by rotating the unit tangent vector along the 98 vessel τ an angle θ tip away from its original direction in an arbitrary plane, as shown in Fig 3. 99 The angle θ tip is chosen from a normal distribution with zero mean and standard deviation σ. 100 An important distinction from previous studies is that the finite extents of the cornea are 101 accounted for: migrating tips are not permitted to leave the simulation domain. A 'soft contact' 102 model is adopted so that tips approaching the boundary of the domain are gradually deflected 103 along the tangent to the bounding surface. Biologically this represents tip cells failing to 104 penetrate the stiffer tissue present on the epithelial and endothelial cornea surfaces, while 'soft' 105 contact is chosen for more robust numerics. The strength of the repulsion from the boundary 106 increases as the boundary is approached according to: where d(x) is the minimum distance to the domain boundaries, d crit is the distance to the 108 boundary at which repulsion is activated and φ max is the dimensionless maximum repulsion strength. The repulsion biases the motion along the tangent to the bounding surface according 110 to: where the· operator produces a random unit tangent. The position of migrating nodes is 112 updated from x(t) to x(t + ∆t) at each time step: where s is the tip velocity and ∆t is the time step size in hours. Fixed values of the repulsion 114 strength φ max = 5 and critical repulsion distance d crit = 25 µm are used for all simulations.  Endothelial tip cells are known to find and merge with other endothelial tips and immature 119 blood vessels during migration, in a process known as anastomosis [23]. There is still 120 uncertainty about the mechanisms by which they meet, but mechanical and chemical guidance 121 are known to contribute [23]. When moving from 2D to 3D models of sprouting angiogenesis it 122 is necessary to define a region within which tips will merge with vessels and each other to allow 123 for the identification of intersections. In this study, a relatively small radius of r ana = 5 µm is 124 used, which is on the order of the vessel radius.

125
VEGF dynamics 126 VEGF dynamics are treated in two different ways in this study. In the first case, the pellet 127 dynamics are ignored and a time independent, spatially varying VEGF concentration field is 128 imposed in the tissue domain. In the second case, the dynamics of VEGF release from a nylon 129 pellet are explicitly modeled.

130
For the first case, a constant VEGF gradient between the limbus and the pellet is applied, 131 given by: where x is a positional coordinate along the geodesic between the limbus and pellet, ε is a small 133 offset from the base, corresponding to the position of the initial vessel, and c p is the VEGF 134 concentration in the pellet, which is assumed to be constant in time in this case. By specifying: 135 the same concentration and concentration gradient magnitudes are maintained at the limbus in 136 all representations. Aside from the Hemisphere, the concentration at the base is 0 nM and at 137 height h + ε (the pellet location) it is c p .

138
For the second case, the transport and decay of the VEGF in the pellet are considered, 139 meaning that c p can now change over time. The dynamics of VEGF in the pellet are not known 140 in detail, but as per Tong and Yuan [7] a high rate of reversible binding to the nylon pellet 141 constituents is assumed. Under this assumption it is possible to derive the following 142 relationship between the total c p and free c f amounts of VEGF in the pellet [5]: where θ ≥ 1 is a dimensionless binding parameter. Free VEGF can decay in the pellet at a rate 144 λ p or leak through the cornea-pellet interface, which has an effective permeability κ p . It is 145 assumed that the VEGF concentration is spatially uniform within the pellet, which has volume 146 Ω p . Balancing mass leads to the following differential equation describing the time rate of 147 change of VEGF in the pellet: where the integral is over the pellet surface ∂Ω and c is the concentration of VEGF in the 149 cornea, at the interface. The initial concentration c p (t = 0) can be determined from the 150 implanted VEGF mass m = 300 ng [5] by: where the VEGF molecular weight MW VEGF is 45 kDa or 45 000 g mol −1 and the factor of 1 1000 152 converts from mol m −3 to M.

153
It is assumed that VEGF diffuses isotropically in the cornea, with a diffusion coefficient 154 D = 2.52 × 10 −7 m 2 h −1 [24,25] and decays naturally at a rate λ = 0.8 h −1 [26]. It is also 155 assumed to enter perfused vessels (and be washed away) and bind to endothelial cells.

156
Combining these processes, we deduce that the dynamics of VEGF in the cornea can be 157 described by the following reaction-diffusion equation: Here κ v = 3 × 10 −4 m h −1 [27,28] is the permeability of vessels to VEGF, R v = 5 µm [29,30] is 159 the assumed vessel radius, c b = 0 M is the amount of VEGF in the blood, assuming it is quickly 160 removed, and ρ and n are respective vessel line and tip densities. The parameter k ec is the rate 161 of VEGF binding per endothelial cell and c 50 is the VEGF concentration at which the rate of 162 binding is half maximal.

163
The extent to which VEGF will pass through the outer cornea layers is not clear, nor 164 whether it will pass through the epithelial layer and into the aqueous humor or the collagen-rich 165 limbus. It is assumed that the rate of such leakage is low, and no flux boundary conditions are 166 imposed on all outer surfaces of the cornea ∂Ω cornea , that is: where n is the inward surface normal. On the cornea-pellet interface the following mass balance 168 is assumed: where c f = cp θ is the amount of free VEGF in the pellet.

170
Eqns (9)-(13) are solved numerically (detailed below), subject to the initial condition of no 171 VEGF in the cornea. In the 2D geometries, the cornea-pellet interface ∂Ω is a line of length w 172 for the planar case or 2πr p for the circle. In the Planar3D geometries it is a rectangle of height 173 T or T p , depending on whether a finite sized pellet is assumed, and width w. In the remaining 174 geometries, the interface is the entire outer surface of the spatially resolved pellet.  Table 1 summarizes the parameter values adopted in this study. Parameter values with sources 177 denoted as 'This Study' are discussed in this section unless previously introduced.

178
A pellet thickness of T p = 40 µm is used in this study, which is less than the T p = 60 µm 179 value reported in Connor et al. [5]. This is to facilitate placement of the pellet in the simulated 180 Hemisphere geometry. The chemotactic sensitivity range χ ∈ [0, 0.5] is chosen to cover extreme 181 cases where the resulting vessel network is not directed towards the pellet and highly directed 182 towards it. The range of the deviation in persistence angle σ ∈ [0, 20] degrees covers cases 183 where straight vessels form, through to cases with vessels with tortuosity similar to that 184 observed in experimental images of the assay. The global time step ∆t = 1 h is chosen to give 185 average segment lengths of 10 µm, which leads to a physically realistic vessel tortuosity. The 186 initial offset from the limbus of ε = 100 µm is in agreement with experimental images [5].

187
The amount of growth factor in implanted pellets is usually known by mass, with a value of 188 300 ng for VEGF reported in Connor et al. [5]. For our pellet volume of 0.0075 mm 3 and 189 VEGF molecular weight of 45 kDa, this corresponds to a pellet concentration of approximately 190 c p = 1330 µM, which is adopted for the dynamic VEGF model.  Simulation results are presented in terms of the 'vessel line density' ρ(x, t), which is defined 206 as the vessel length per unit volume, and 'tip density' n(x, t), which is the number of migrating 207 tips per unit volume. Densities are calculated on structured grids, with values averaged over 208 grid cells that are equidistant from the limbus. To reduce noise caused by the sampling of 209 discrete vessels and tips onto the grids, two Gaussian smoothing passes are applied to ρ and n 210 before further processing [33].  leading to a reduced number of tips and greater confinement of tips toward the advancing front. 219 In 3D, the presence of multiple vessels through the cornea thickness is evident. In the Circle2D 220 case, there is a tendency for tips to move together as the center is approached, this focusing 221 effect being due to the domain geometry. concentration is localized to a line of length w on the interface. The Planar3D case, with finite 243 pellet width, has a noticeably lower VEGF concentration than the 2D case, due to the pellet 244 thickness T p being smaller than that of the cornea T . In the circular domains, the situation is 245 reversed, with higher VEGF concentrations in the 3D domain due to a greater surface density 246 at the cornea-pellet interface. A higher concentration is observed in the Hemisphere for the 247 same reason.

248
Over time the VEGF in the pellet depletes, with the decay term in Eqn (9)  comparable with the Hemisphere, due to focusing of the vascularized region. When the pellet is 258 moved closer to the limbus the general trend is for an increase in the maximum tip density. In 259 relative terms, the tip density increases most in the Hemisphere and Planar3D geometry with 260 finite pellet, both by a factor of 1.7, although in absolute terms the density in the 261 Planar3DFinite geometry is approximately half that of the Hemisphere. In contrast, in the 262 circles, the tip density is reduced by a factor of 1.25 as the pellet is moved towards the limbus. 263 This is a geometric effect, due to the breaking of symmetry as the pellet is moved away from 264 the center of the circle. The dynamics of the maximum densities are similar to those shown in 265 Fig 4C, with the rate of density increase being greatest in the circular domains.  The location of the maximum tip density and half-maximum tip density are useful additional 274 metrics in cases where they can be measured. As shown in Fig 6C, these metrics are more can be difficult to measure experimentally. This is because endothelial tip cells are not obvious 289 at typical imaging resolution, while line density measurements are subject to potential errors as 290 multiple vessels may overlap through the cornea thickness [7]. In contrast, it is easier to 291 estimate the 'vascularized fraction' or volume of the domain with vessels divided by total 292 domain volume directly from images. In the present study this metric is calculated by 293 accumulating the volume of the cells in the structured grid used in the calculation of densities 294 that are occupied by vessels and dividing by the volume of all cells in the grid.

301
Summarizing the results in Figs 4, 5, 6 and 7, it is predicted that the small, but finite, thickness 302 of the cornea can have an important effect on vessel network formation relative to 2D assays 303 and models, reducing the likelihood of vessels anastomosing. The positioning of the pellet 304 affects the resulting neovascularization in the Hemisphere, with increased focusing observed as 305 the pellet is moved towards the limbus. Circular domains lead to increased anastomosis relative 306 to planar domains and the hemispherical cornea. Circular domains also have a markedly 307 different sensitivity to pellet positioning.

308
The curvature of the cornea is predicted not to affect neovascularization, with good 309 agreement between the Planar3D domain with a finite pellet and Hemisphere simulations.  The 'location of the vascular front' is predicted to be robust to the choice of domain, 319 suggesting that it could be used to relate in silico predictions or in vitro observations to in vivo. 320 The 'vascularized fraction' is predicted to be sensitive to domain choice, in particular the 321 manner in which VEGF is introduced to the system. This may help to relate observations in 322 planar domains to those in the in vivo assay.

323
In the present work, a simplified model of neovascularization is adopted, with the primary 324 focus being on comparing predicted neovascularization across different geometries, which has 325 not been performed for the cornea before. By varying phenomenological parameters, such as regression, which may reduce overall vessel densities and average vessel lengths and ultimately 331 halt the progression of the moving front [40], ii) a more detailed model of endothelial cell 332 proliferation and migration [6], which would remove the need for a phenomenological migration 333 speed s and may lead to a front velocity that is not constant in time, iii) models of tip 334 attraction [40] and mechanical interaction with extracellular matrix [6], which may encourage 335 extra anastomosis by guiding tips to locate each other and iv) modelING detailed feedback 336 between a metabolically active tissue and the vasculature [40]. 337 Further, there is scope for more detailed parameter studies, such as an investigation of the 338 relationship between pellet location and 'vascularized fraction' and the sensitivity of predictions 339 to the anastomosis radius r ana . More quantitative comparison with experimental measurements 340 of the kind in Tong and Yuan [3] and Connor et al. [5] is now possible, due to the focus on 341 modeling the actual geometry of the cornea-pellet system in the present study. Such In this study we use a newly developed 3D, off-lattice mathematical model of neovascularization 351 in the cornea to study how: i) the geometry of the cornea-pellet system in the micropocket 352 assay affects vessel network formation and ii) which metrics of neovascularization are most 353 sensitive to geometrical differences between typical in silico, in vitro and in vivo tissue domains. 354 We predict that:

355
• 2D and circular domains lead to increased anastomosis, even relative to the thin cornea 356 geometry,

357
• predictions of neovascularization are highly sensitive to the geometrical treatment of the 358 VEGF-containing pellet,

359
• measuring the distance of the growing vascularized front to the limbus leads to 360 predictions that are robust to differences in domain choice,

361
• vascularized fractions can serve as a useful proxy for densities of migrating tips or vessel 362 line densities, which are more difficult to measure. These metrics can better distinguish 363 underlying vascular patterning than the distance to the limbus,