In vitro and computational modelling of drug delivery across the outer blood–retinal barrier

The ability to produce rapid, cost-effective and human-relevant data has the potential to accelerate the development of new drug delivery systems. Intraocular drug delivery is an area undergoing rapid expansion, due to the increase in sight-threatening diseases linked to increasing age and lifestyle factors. The outer blood–retinal barrier (OBRB) is important in this area of drug delivery, as it separates the eye from the systemic blood flow. This study reports the development of complementary in vitro and in silico models to study drug transport from silicone oil across the OBRB. Monolayer cultures of a human retinal pigmented epithelium cell line, ARPE-19, were added to chambers and exposed to a controlled flow to simulate drug clearance across the OBRB. Movement of dextran molecules and release of ibuprofen from silicone oil in this model were measured. Corresponding simulations were developed using COMSOL Multiphysics computational fluid dynamics software and validated using independent in vitro datasets. Computational simulations were able to predict dextran movement and ibuprofen release, with all of the features of the experimental release profiles being observed in the simulated data. Simulated values for peak concentrations of permeated dextran and ibuprofen released from silicone oil were within 18% of the in vitro results. This model could be used as a predictive tool for drug transport across this important tissue.


Introduction 1
The ability to produce models that mimic tissue biology and physiology is crucial in the development of drug 2 delivery systems. Currently, much of the pre-clinical work is conducted in animals, the limitations of which 3 are well-reported and which may be the cause of the 86% of drugs that fail in clinical testing (1). The use of 4 in vitro models means that human tissue can be used, minimising species differences as well as ethical 5 concerns. Although lacking the complex environment of a living body, they can be used to generate reliable 6 data on drug transport and toxicology (2-4). Computer science and simulation of pharmacokinetic behaviour 7 is becoming an integral part of pharmaceutical research and development due to the significant reductions in 8 cost and time they can offer. In silico tools can, in conjunction with complementary empirical measurements, 9 contribute to the optimisation of novel drug delivery systems, and ever-increasing computing power has 10 allowed more sophisticated models to be developed. The combination of accurate, validated in vitro and in 11 silico models has the potential to revolutionise the development of drug delivery technologies by providing 12 rapid, reproducible and human-relevant data (5). 13 Epithelial barriers are particularly important in drug delivery, as one of the primary functions of epithelium is 14 to act as a barrier, and certain drugs must cross the barrier to reach the target tissue. Tissues such as skin (6), 15 intestinal (7) and pulmonary epithelium (8) have all been studied in the development of drugs and chemicals, 16 and the complexity of the models of these barriers has improved dramatically in the last few decades, from 17 simple two-dimensional, multi-well plate cultures to intricate microfluidic culture chips. Many epithelial 18 barrier models, such as those for skin, are based on an air-liquid interface,. In the case of fluid-fluid 19 interfaces, many are designed as static culture systems. These static culture systems allow the determination 20 of parameters such as barrier permeability and diffusion coefficients, which are of crucial importance when 21 studying the movement of potential therapeutic treatments across epithelial tissues. They do not, however, 22 provide realistic, time-dependant development of concentration gradients across the model as they do not 23 mimic the many dynamic factors associated with these physiological barriers. With the rapid rise in 24 microfluidic technologies for cell culture that has occurred over recent years, the creation of dynamic in vitro 25 models has become more accessible (9). 26 There are many eye diseases that require treatment with pharmacological agents. In contrast to the front of 27 the eye, where drug delivery can often be achieved by topical application of eye drops and ointments, many 28 eye diseases of the posterior segment of the eye require the delivery of drugs directly into the vitreous. The 29 drugs may be required to treat acute infection or inflammation, or to treat a chronic condition such as age-30 related macular degeneration or diabetic retinopathy. In the latter case, repeated intravitreal injection or use 31 of implantable drug delivery devices are the preferred methods to achieve therapeutic levels of drug over the 32 extended periods required. Due to the cost and invasive nature associated with repeated intravitreal injection, 33 as well as the potential for sight-threatening complications, much effort has been directed towards 34 developing implants that can deliver drugs over extended periods (10). We have recently developed 35 technology to achieve extended release of drugs from silicone oil tamponades (11,12). Silicone oil and gas 36 tamponades are used to replace the native vitreous humour in the treatment of sight-threatening retinal detachments. They inhibit the flow of aqueous fluids into the subretinal space, exclude inflammatory factors, 38 and support the retina as tears heal (13,14). Silicone oils are the only medical devices licenced for long-term 39 use as tamponades. Drug release from silicone oil tamponades would release pharmacological adjuncts to the 40 surgical treatment with the aim of reduce complications caused by scarring conditions such as proliferative 41 vitreoretinopathy and proliferative diabetic retinopathy. The ability to predict the release of drugs from such 42 devices, as well as understand factors that influence clearance from the posterior cavity, is crucial if over-43 and under-dosing is to be avoided. 44 Drugs administered intravitreally, whether via injection or an implant, are cleared via two routes, either 45 anteriorly or posteriorly ( Figure 1A). Nearly all compounds can be eliminated via the anterior route. 46 Anteriorly, once the drug has diffused across the vitreous, it enters the aqueous of the posterior chamber 47 where it is transported into the anterior chamber and cleared, either through the Schlemm's canal or into the 48 uveal blood flow (15). Alternatively, posterior elimination occurs through permeation across the posterior 49 blood-ocular barriers such as the outer and inner blood-retinal-barrier (BRB). This route requires either 50 adequate passive permeability of agents, which is generally only applicable to very low molecular weight or 51 lipophilic substances, or active transport of molecules by the cells present in the barrier. For this reason, 52 larger or hydrophilic molecules generally have longer half-lives in the vitreous as they are not able to move 53 as freely through the vitreous or pass across the retina (16). 54 In vitro models have often been used to investigate the permeability of drugs and compounds across the 55 retinal pigmented epithelium (RPE) or to model diseases of the outer blood-retinal barrier (OBRB). In order 56 to obtain a physically accurate model of the OBRB, the main anatomical structures (retinal pigment 57 epithelium, Bruch's membrane and choroid), as well as physiological conditions, such as flow, need to be 58 incorporated. Although many of the in vitro models of the OBRB that have been reported include the main 59 anatomical structures within the tissue (17-20), few include a flow mechanism to model the blood flow 60 within the choroid and to avoid the formation of an unstirred water layer (21). This flow mechanism is 61 particularly important when investigating the clearance of drugs across the BRB, as it causes systemic by, for example, incorporating flow, leaves scope for this to be investigated further.
At time intervals of 1, 3, 8 and 24hr, 50µl samples were taken from the receptor compartment and replaced 143 with fresh PBS solution. Samples were stored at 4⁰C and protected from light until the time points had been 144 completed, and the fluorescence intensity of the samples was read in a microplate reader at λex (excitation 145 wavelength): 485nm, λem (emission wavelength): 535nm. Flux was determined from the slope of the linear 146 portion of the curve. Calibration curves (R 2 ≥0.997) for each molecule were made using a serial dilution of 147 50µg/ml dextran solution. Average measurements for blanks (0µg/ml of dextran) were subtracted from the 148 standards and the unknown samples. Concentrations of the unknown samples were determined from the 149 calibration curves and each sample was repeated in triplicate (n=3). 150

Measurement of drug release 151
An experimental concentration of 1mg of ibuprofen (ibu) in 1ml of 1000c.st silicone oil (Fluoron GmbH) 152 was used. The ibu-SiO was stirred for 72 hours in a sealed flask and then filtered in a Class II biological 153 safety hood to sterilise. In a 24-well plate, 1mL of ibu-SiO was syringed on top of 500µl of PBS. At defined 154 time points, 100µl samples were taken from the PBS, transferred to UV transparent cuvettes (UVette, 155 Eppendorf) and the time-dependent increase in ibuprofen concentration measured by UV-Vis spectroscopy 156 (n=3). 157

Determination of Diffusion coefficients 158
The diffusion coefficients (D) were calculated using the Stokes-Einstein equation, using the appropriate 159 solvent viscosity for either DMEM-F12, PBS, water or SiO (µ), the apparent radius (r) of either the dextran 160 or ibuprofen molecule (36), temperature (T) and Boltzmann constant (k): 161 D=kT/6πμr 162

Drug clearance study 163
The Kirkstall QV600 system was arranged in a single-pass series fluid circuit using a peristaltic pump to 164 control the fluid flow inlet rate ( Figure 1C). Either acellular or cell-seeded cell culture inserts containing 165 400µl of tracer solution were inserted into the QV600 chamber to allow flow across the receptor side of the 166 membrane. The cell culture inserts have a smaller diameter than the QV600 chamber. In order to create a 167 seal between the donor and receptor compartments, a silicone O-ring was placed around the exterior wall of 168 the insert before it was placed in the chamber. The system was incubated at 37⁰C and phosphate buffered 169 saline solution was perfused through the receptor chamber at a constant flow rate for 24 hours (8 hours for 170 2ml/min experiments). The flow rates used were 20µl/min, 200µl/min, 400µl/min. 171 FITC-conjugated dextran (4kDa) was dissolved in DMEM-F-12 supplemented with 10% foetal calf serum to 172 50µg/ml in the donor chamber. The systems were set up and perfused for 1hr, after which the flow was 173 stopped and the tubes were clamped. The solution in the receptor chamber was completely removed and 174 homogenised, and 50μl samples were taken. The systems were cleaned and reset, and the protocol was 175 repeated for increasing periods of time (1hr time increments). The fluorescence intensity of the samples was read in a microplate reader at λex: 485nm, λem: 535nm. Concentration of the unknown samples was 177 determined from the calibration curves and each sample was repeated in triplicate (n=3). 178

Measurement of ibuprofen release from silicone oil under flow conditions 179
The QV600 chamber and peristaltic pump were assembled as described previously. The system was primed 180 with 30mL of sterile PBS. 2mL of 1mg/mL ibu-SiO was added directly on top of the PBS. The system was 181 incubated at 37⁰C and PBS was perfused through the receptor chamber at a constant flow rate for 24 hours (8 182 hours for 2ml/min experiments). The flow rates used were 20µl/min, 200µl/min and 2ml/min. At set time 183 intervals, the ports were clamped shut and the volume of PBS beneath the ibu-SiO removed using a 25-gauge 184 needle. A 25G needle allows the PBS to be removed but the viscosity of the oil prevents its withdrawal 185 through the needle. The solution was homogenised and 50µL samples taken. UV-Vis was used to determine 186 the ibuprofen concentration in the samples as previously described. 187

Computer Model 188
Geometry of the QV600 chamber and grid generation 189 were dependent on whether a representative cell monolayer was included. 196 The mesh was generated using the commercial software, COMSOL Multiphysics. For the single phase, 197 dextran transport studies, the mesh comprised free triangular elements with boundary layers at the no slip 198 walls. The total number of elements in the mesh was 38267. For the ibuprofen release studies the mesh 199 consisted of 2258 free triangular elements in a moving mesh system to model the flow of two immiscible 200

Governing equations 202
The incompressible Navier-Stokes equations along with a species transport equation were solved in order to 203 obtain the velocity and concentration fields across the models. 204 ci.. The diffusion coefficient, D, which was determined from the in vitro studies, is specific to each species.

217
These are presented in Table 1. As the drugs used in this study are not reacting with the cells and the cell 218 layer acts only as a barrier, the right-hand side of the equation is zero. 219 Two models were investigated in this study: one to study the passage of dextran molecules across the 220 NH3_ePTFE_M membrane and the other to study the release of ibuprofen from silicone oil. These two 221 models are described separately below. 222

Dextran Transport Model 223
The model itself consists of two parts: a laminar flow interface to compute the velocity flow and pressure 224 fields of the single-phase fluid flow, and a transport of diluted species interface. COMSOL provides this 225 interface to calculate the concentration field of a dilute solute in a solvent, i.e. the fluorescently labelled 226 dextran solutions diluted in culture medium. The model was run to simulate both the absence and presence of 227 cells using alterations in both the geometry and permeability boundary conditions that are described below. 228 An additional physics node was included to model the transport of FD through the membrane into basolateral 229 medium flow. This model accounts for the dissipation of kinetic energy experienced by the fluid moving 230 through a porous matrix through means of viscous shear. In COMSOL, this is implemented using the fluid 231 and matrix properties node which uses the Brinkman equations: 232 where κ is the permeability of the membrane to the fluid and εp is the porosity (Table 1). This node was only 234 applicable to the membrane part of the model and were therefore only applied to that domain. 235

Boundary Conditions 236
A no slip boundary condition was applied to the walls of the geometry: 237 = 0 238 The inlet was applied to the left-hand wall of the inlet tube. The velocity field, v, for the inlet was defined 239 where Q is the flow rate and A is the cross-sectional area of the inlet tube. A range of inlet flow rates was 242 investigated to coincide with the flow rates used in the in vitro experiments: 20µL/min, 200µL/min, 243 400µL/min, and 2mL/min. The outlet was applied to the right-hand wall of the outlet tube. The outlet 244 condition was a zero pressure (p) condition: 245 Two transport of diluted species nodes were used in the transport of dextran simulations: one for the donor 247 and receptor domains, and one for the membrane domain. This separate node for the membrane domain 248 allowed the difference in diffusion in that domain to be accounted for. The second transport of diluted species node applied similar boundary conditions as above but across the 259 membrane domain, therefore accounting for the difference in diffusion. A no flux condition was applied on 260 the exterior wall boundaries of the membrane. The same function for the pointwise constraint was applied to 261 the boundaries that were shared with the other two domains. The dimensions of the membrane domain were 262 altered depending on whether the presence of the cells was being modelled or not. In the presence of cells, 263 the membrane domain was 70µm in height and the appropriate diffusion coefficient of the domain was used, 264 as described in Table 1. In the absence of the cells, the domain was reduced to 50µm and the diffusion 265 coefficient was altered to account for their absence. 266

Ibuprofen Release Model 267
The second model explored ibuprofen release from SiO. With the differences that the ibuprofen release 268 studies required, the model was altered to simulate the interaction of two immiscible fluid phases: the 269 aqueous PBS phase and the 1000 c.st silicone oil phase. This model also removes the membrane domain as 270 the movement of drug was directly from the oil into PBS. 271 To model the two-phase nature of the model, an additional moving mesh mechanism was used. The laminar 272 flow moving mesh physics node in COMSOL solves the same equations for velocity and pressure fields, but 273 also tracks the movement of the interface between two immiscible fluids by allowing deformation of the 274 mesh during the solution. A free mesh deformation was prescribed to the domains either side of the interface. 275 The inlet and outlet tubes of the geometry were prescribed a fixed mesh as only one of the fluid phases 276 moved through these regions. The mesh is also prescribed zero displacement at the exterior boundary walls 277 to prevent collapsing of the solid wall boundaries. The walls that were in contact with the fluid-fluid 278 interface were prescribed free deformation of the mesh parallel to the exterior wall boundaries, but with zero 279 perpendicular displacement, again to prevent collapse of the solid exterior walls. 280 The additional equation solved for in the laminar flow moving mesh interface was the Navier slip equation. 281 This condition was applied to the boundaries that were in contact with the fluid-fluid interface, and is 282 appropriate for the two-phase flow model. This condition adds a frictional force, Ffr, at a stationary wall, 283 which allows the interface to move against the wall. 284

= − 285
where β is the slip length, which was a function of the element size of the mesh, μ is viscosity and u is the 286 velocity vector. The fluid-fluid interface node also takes into account the interfacial tension of the two fluids, 287 σ = 50mN/m (37), and the contact angle between the wall and the fluids, θw = 1.3 rad (38). 288 Due to the nature of the two-phase model, a stationary solution for the velocity field could not be solved 289 because of the movement at the interface, therefore a time-dependent solution was obtained over 9 seconds, 290 at which point the flow stabilised. This stabilised flow was used as the velocity field input for the transport of 291 diluted species solutions. 292 The solution for the transport is simpler for the release of ibuprofen because the concentration species only 293 moves through two domains and not the fluid matrix domain, therefore the movement is purely diffusion and 294 convection in the two different fluids. An additional expression for the partition coefficient (Pp) is included 295 in the mass transport between the two phases. 296 where Jd is flux, Ps is permeability and cd and cr represent the concentration in the donor and receptor domain 298 respectively. 299 300

Boundary Conditions 301
A no-slip condition was applied to the exterior boundary walls. The walls in contact with the fluid-fluid 302 interface were assigned a Navier slip condition. This condition allows the fluid-fluid interface to move along 303 the wall. Additionally, the top boundary of the oil phase, parallel to the fluid-fluid interface, was assigned a 304 slip condition. This slip condition allows the deformation of the mesh to continue throughout the phase 305 whilst still applying a no penetration condition, meaning the model allows the movement of the mesh 306 without fluid leaving that domain. 307 The inlet and outlet conditions were as described previously and the same flow rates were investigated as 308 with the dextran transport studies (20µL/min, 200µL/min, 400µL/min and 2mL/min). A volume force was 309 also implemented across the entire geometry to account for gravity in the system. 310 Two transport of diluted species nodes were used to investigate the release of ibuprofen from the silicone oil: 311 one for the oil phase and one for the aqueous phase. The appropriate diffusion coefficients (Table 1) were 312 applied to each fluid domain and a pointwise constraint was applied at the fluid-fluid interface. This 313 pointwise constraint takes into account the concentration at the interfaces and solves for the mass flux across 314 that boundary using a function of the concentration gradient and the partition coefficient. The accepted error 315 between the computer models and experimental data was set by two boundaries: <10% was considered to be 316 good agreement and <20% was considered to be acceptable agreement. 317

Material properties 318
1000 c.st silicone oil is often used by surgeons because its viscosity makes it easily injectable, a motivation 319 for its use in this study. SiO has a lower density than water, and therefore floats on it, but this oil has a higher 320 viscosity. These properties are shown in Table 2. Rheological evaluation of water, PBS and silicone oil was 321 done using a Rheosense Inc µVISC rheometer (Rheosense Inc., USA) with a 100N load cell. The results 322 showed negligible differences between water and PBS in terms of viscosity and density, therefore water was 323 used as the aqueous fluid of interest in the computer model. A handheld density meter (Anton Parr) was used 324 to measure the density of the fluids at 37°C. 325

Grid Independence Studies 327
To demonstrate grid independence, simulations were run with varying degrees of mesh refinement (dextran 328 model between 18,000 and 108,814 elements, ibuprofen model between 300 and 4000 elements). The 329 velocity at two points in the geometry was measured (Figure 2A, 2E) and compared as the mesh was refined. 330 The number of elements used in each model was justified by using the mesh refined to within 5% agreement 331 between both points and the highest resolved mesh ( Figure 2G, H). This was chosen as an acceptable error 332 limit that reduced the computational time whilst maintaining sufficient accuracy. For the dextran model, this 333 was 38267 elements. For the ibuprofen model, this was 2258 elements.  (Table 1) were comparable with those seen in a similar experiment by Mannermaa et al. (19). In that study, 359 the authors compared the transport of drugs through a static, ARPE-19-based model and bovine RPE tissue, 360 finding similar transport trends. The simulation shows the concentration gradient of dextran in the donor 361 compartment to decrease with increasing flow rate after 24 hours. As expected, dextran is less readily cleared 362 across the barrier when the presence of cells is included in the simulation (Figure 4). The results of the model 363 were then validated using data from complementary in vitro experiments. The numerical model shows 364 agreement with the in vitro data in both the acellular and cell seeded experiments (Figure 5). At the higher 365 flow rates, the simulated results are able to mirror the change in release exhibited in vitro, which shows a 366 shift to a burst release response followed by an exponential decay in concentration over time, and there 367 appears to be no correlation between flow rate and simulation accuracy (42). The simulation of dextran 368 transport was able to predict the maximum concentration (Cmax) observed in that chamber to within 5% of the 369 acellular experimental data. The introduction of cells to the system increased the error observed in Cmax but 370 still to within 18% of the experimental data. Similar studies, which simulated permeability of different 371 molecular weight FITC-dextrans in a static set-up and across collagen or agarose gel, showed an increase in 372 error (between approximately 4% and 46%) between their simulated and experimental results with increasing molecular weight (43). Here, we have only presented data for the transport of 4kDa FITC-dextran; other 374 sizes were also investigated (40kDa and 70kDa FITC-dextran) with 4kDa and 40kDa producing similar 375 errors, but 70kDa showed increased error in comparison. The simulation does not take into consideration the 376 biological effects of culturing cells under flow might have on the barrier functionality of the ARPE-19 cell 377 monolayer. There are studies which have investigated biological effects on epithelial tissues in computer 378 simulations, for example modelling inflammatory effects on intestinal epithelium in necrotising enterocolitis 379 (44), and investigating links between epithelial morphogenesis and cancer mutations (45). A combination of 380 computational fluid dynamics modelling such as in this study and a more computational biology approach to 381 investigate cell dependent changes in transport and clearance of molecules could further improve the 382 agreement between the experimental and simulated data in cell seeded simulations. 383

Ibuprofen Release Studies 384
Tracking the fluid-fluid interface is important when considering the concentration distribution of drugs 385 across the two fluids. An adaptive mesh was used to simulate the interaction between the silicone oil phase 386 and the aqueous (PBS) phase. Based on the interfacial tension and wall contact angle of the two fluids, the 387 simulation shows that a meniscus is formed between the two phases over time until a steady state is reached 388 at approximately 1.4 seconds ( Figure 6). 389 In the two phase system (oil and aqueous), the simulation shows two distinct flow fields within each phase 390 ( Figure 7). The flow field in the oil phase, however, does appear to be influenced by the flow rate of the 391 aqueous phase. The flow fields formed at 20µL/min, 200µL/min and 400µL/min show similarities to those 392 formed in the single phase, membrane system. At 20µL/min, however, the low flow rate inlet stream within 393 the aqueous phase appears to bounce off the oil phase and creates a ripple within the primary flow stream. 394 The interaction between the main flow stream and the oil phase also creates two recirculating streams within 395 the oil domain itself. This phenomenon is observed for each of the flow rates studied, with the split in the 396 two streams shifting towards the inlet as the inlet flow rate increased. Unlike the 2mL/min flow field in the 397 To model the release of ibuprofen from SiO, the physics controls of the computer model was redesigned to 405 allow the interaction of two immiscible fluid phases within the QV600 chamber. For this ibu-SiO model, the 406 moving mesh method was applied. COMSOL Multiphysics reports that this method provides the best results 407 when tracking the interface between the phases is important, and also allows mass transport across the 408 interface, which is difficult to implement using other methods (46). At 20µL/min, the flow entering the main chamber bounces of the bottom of the meniscus formed by the oil phase and creates a ripple in the main

Author Contributions 480
AED carried out lab work, computational modelling, data analysis and drafted the manuscript. RLW 481 participated in the design of the study and critically revised the manuscript. GL critically revised the 482 manuscript. SRP participated in computer model design, data analysis and critically revised manuscript. 483 VRK conceived the study, designed the study, coordinated the study and helped draft manuscript. All authors 484 gave final approval for publication and agree to be accountable for the work performed therein. 485

Data Accessibility Statement 486
The datasets supporting this article are available online (42). DOI: 10.17638/datacat.liverpool.ac.uk/908. 487 Input geometry, red dots indicate points 1 and 2 for grid independence study, and (F) mesh generated for 629

Figure Legends
Kirkstall QV600 chamber used in ibuprofen release simulations. This mesh comprised 2,258 free triangular. 630 (G) Percentage difference in velocity at two points in the centre of the chamber compared with the finest 631 mesh for dextran transport (finest mesh: 108,814 elements). (H) Percentage difference in velocity at two 632 points in the centre of the chamber compared with the finest mesh ibuprofen release (finest mesh: 3481 633 elements). An acceptable mesh density was deemed to produce <5% error. 634   Experimental data presented as mean concentration ± 1 SD, n = 3. There is good agreement with respect to 647 the concentration at all flow rates, including the change of behaviours from burst release to exponential 648 decay, with differences within 5% for acellular and 18% for cellular experiments. 649