Multiscale modelling of drug transport and metabolism in liver spheroids

In early preclinical drug development, potential candidates are tested in the laboratory using isolated cells. These in vitro experiments traditionally involve cells cultured in a two-dimensional monolayer environment. However, cells cultured in three-dimensional spheroid systems have been shown to more closely resemble the functionality and morphology of cells in vivo. While the increasing usage of hepatic spheroid cultures allows for more relevant experimentation in a more realistic biological environment, the underlying physical processes of drug transport, uptake and metabolism contributing to the spatial distribution of drugs in these spheroids remain poorly understood. The development of a multiscale mathematical modelling framework describing the spatio-temporal dynamics of drugs in multicellular environments enables mechanistic insight into the behaviour of these systems. Here, our analysis of cell membrane permeation and porosity throughout the spheroid reveals the impact of these properties on drug penetration, with maximal disparity between zonal metabolism rates occurring for drugs of intermediate lipophilicity. Our research shows how mathematical models can be used to simulate the activity and transport of drugs in hepatic spheroids and in principle any organoid, with the ultimate aim of better informing experimentalists on how to regulate dosing and culture conditions to more effectively optimize drug delivery.


Microscale Transport
The mathematical representation of microscale drug transport across a cell membrane can be studied with a simple model that considers the processes governing drug concentration dynamics in two phases, inside and outside the cell, with a permeable barrier in-between. Assume diffusion occurs at different rates inside ( " ) and outside ( # ) of a spherical cell of radius and that the drug is metabolised within the cell. The drug concentration ( ) dynamics inside the cell are given by the partial differential equation (PDE): with scaling where 012 is the maximum metabolism rate and 0 represents the drug concentration at which metabolism is half maximal. Note that the metabolism of the drug is assumed to occur with Michaelis-Menten kinetics, which is relevant for enzyme-mediated biochemical reactions but here prohibits the derivation of an analytical solution at the steady state. Assume that outside the cell drug transport is governed by diffusion processes only: where = # / " . For simplicity, assume that the problem is radially symmetric, drop the tildes and convert to spherical coordinates for a 1D representative model with respect to the radius, .
We impose the following boundary conditions for the cell centre ( = 0) and a distance away from the cell = 012 : where B CDE is a constant supply term to be prescribed. Assume that the flux at the sphere boundary is equal such that mass is conserved, i.e., where " and # are used to distinguish between interior and exterior drug concentrations at the cell membrane boundary. A further boundary condition must be specified at the cell membrane boundary in order to solve the coupled PDE system and investigate the effects of different means of drug transport.

Passive diffusion
To determine the boundary condition describing the flux of drug into the cell due to passive diffusion, consider an additional compartment, i.e., the cell membrane. It is assumed that within this compartment, drug transport is determined solely by aqueous diffusion. Since the thickness of the membrane (~5-10 nm) is much smaller than the cell radius (~10-20 µm) and surrounding cell space, it is assumed that the drug diffuses across the space between the lipid barriers of the membrane relatively quickly compared to transport outside. Therefore, we assume that there is a valid quasi-steady state assumption to be made at either side of the membrane such that drug concentration can be assumed to be constant at the lipid barriers on this quick timescale. Mathematically, we can represent this as a thin membrane compartment of width in which the drug concentration ( G ) is transported across the space via diffusion (at a rate G ) with Dirichlet boundary conditions: We can solve equation (S9) via integration to give where and are constants that are determined using the boundary conditions (S10)-(S11) such that Therefore, by substitution into (S12), we acquire the solution for ≪ 1. The inward flux at the cell membrane can then be found by Fick's law: where G / = represents the permeability coefficient, proportional to the rate of intra-membrane aqueous diffusion and inversely proportional to the thickness of the membrane. Since ≪ 1 we can show that the outer membrane flux can also be derived such that and the inner and outer fluxes are equal as → 0 and thus, " " = ( # − " ) = # # , = 1 . (S18)

Numerical solution
We can solve the system numerically using the method of lines and gears whereby the following finite difference approximations in the spatial dimension reduce our PDE problem to an ODE problem. We apply central difference formulae for 1 st and 2 nd order spatial derivatives: for = 1. . ( = 0. .1) where we use notation \]^= ( + ∆ , ) such that an increase by 1 in subscript corresponds to a radial increment in ∆ as defined by the discretisation of the mesh (spatial domain). Similarly our PDE model for exterior dynamics is reduced to for = + 1. . + 1 ( = 1 + ∆ . . 012 ). We now inspect our boundary conditions in order to determine special case boundary values.

Sphere centre boundary
In the case where = 1 we have and therefore a singularity at = 0. In the limit, by use of the Neumann boundary condition in equation (S6) and L'Hôpitals rule. Consequently, Thus we need to determine the value of the node at = 0, i.e. = 0 − ∆ . In order to do this we apply the Neumann boundary condition in equation (S6). Therefore,

Phase interface boundary
For the interface boundary of the sphere we have a discontinuity and the following equations for the concentrations either side of the boundary: where special boundary values have been highlighted with accents for extra clarity (note in passive diffusion uninhibited by any membrane permeation, k j]^= k j ).
When = , we define \]^= k j]^ using the equal flux boundary condition (see equation (S8)) and one-sided finite difference approximations and consequently we have To determine k j we use the transport boundary condition in equation (S18): Therefore, following substitution into equation (S34), where / = / " and substituting equation (S38) into equation (S37),

External boundary
and we determine the value of at = + 1, i.e. = 012 . In order to do this, we apply the Dirichlet boundary condition in equation (S7): Therefore,

Model simulations
The systems of ODEs derived in section 1.1.1 were numerically integrated using MATLAB R2017b software to illustrate typical solutions (parameters provided in the legend for Figure 1)  where ~Bt{ ( , ) is the model solution at the original/default resolution (1001x1001 mesh) and •Bt{ ( , ) is the model solution at an increased mesh resolution where the spatial discretisation is increased 10-fold).

Carrier-mediated transport
We here reiterate a brief derivation of the simple carrier model used to derive the membrane flux boundary condition for clarity. It is assumed that the drug substrate (external, # , or internal, " ) can reversibly bind to the transporter (facing the exterior, , or interior, ′, of the cell), with first order, mass-action kinetics. The bound transporter complexes are given by the variables and ƒ for outward and inward facing transporters respectively. It is also assumed that the transporter (bound or unbound) undergoes a conformational change with first order kinetics to change position such that the binding site is facing either the exterior or interior of the cell (see schematic in Fig1A). The above processes are described in the system of ordinary differential equations below: where ˆ terms represent conformational changes, ‡ terms represent binding/unbinding and the amount of receptor is conserved, i.e., ′ + + ′ + = b (constant). The drug substrate flux is given by the difference between interior dissociation of bound substrate and the association of unbound substrate with inward facing receptors (and is equal to the difference between the association of unbound substrate to outward facing receptors and substrate dissociation at the cell exterior). Thus the flux is, defined to be positive from outside to inside. It is assumed that the processes of binding and conformational changes are fast relative to the spatiotemporal drug concentration dynamics of the model at the cellular scale. Therefore we can find the steady state flux by setting the left-hand sides of equations (S43)-(S46) equal to 0 and solving to derive the 4 state variables in terms of # , " and the rate constants, subject to total receptor concentration, b . By substitution into equation (S47) we acquire the steady state flux: .

Numerical solution
The carrier-mediated transport model was solved numerically using finite difference approximations as before. However, the new boundary condition in equation (S54) required the introduction of modified boundary values in the ODE approximations at the phase interface boundary, With terminology analogous to the passive diffusion model, we use one-sided finite difference approximations and the carrier-mediated flux boundary condition in equation (S54) to determine k j : where / b = b / " . Therefore, following substitution into the equal flux equation (S34), we can derive a quadratic equation for k j]^, and numerically integrate the full system of ODEs as described in section 0.

Diffusion of small molecule drugs
From a sample data base of 321 drugs, we identified several important physicochemical properties including molecular weight and density (Kyffin, 2018). These values allowed us to formulate an estimated range of likely diffusion coefficients for a wide range of drugs with a physically accurate range of weights and densities by employing the Stokes-Einstein equation describing the diffusion of spherical particles through a liquid with low Reynolds number, where ‹ is Boltzmann's constant, represents temperature (assumed to be the physiological value, 310.15 K), is viscosity (assigned as 6.913×10 -4 kg m -1 s -1 , the dynamic viscosity of water at 310.15 K) and is the particle radius (m). To calculate the radius, we assume drugs can be represented spherically and use MW and density ( ) data: By implementing the formulae in (S61) and (S62) for the drug data base (MW ~100-1,200; density ~0.6-2.6 g/m 3 ), we were able to identify a feasible diffusion coefficient range of approximately 5×10 -10 to 1×10 -9 m 2 /s. At early times, we assume that passive diffusion can be represented by the following system where the rate constant \u represents the transport of drug into the cell:

Permeability as a function of lipophilicity
where ˆtšš is the amount of drug (units of moles) in the cellular compartment (expressed per 10 6 cells) and b is the substrate or media concentration (moles dissolved in 0tv = 400 µL of media in the Menochet et al. study), which we assume to be an approximate constant external supply at early times, i.e., equivalent to dose concentration. The passive diffusion clearance, v\•• , is defined as the slope of the uptake rate against concentration, for media substrate ¡¢ and dimensions, where [V] = volume units. By comparison with equation (S64), at early times, we have where v\•• has units of µL/min/10 6 cells in the Menochet et al. study.
In order to translate this uptake-related parameter into our spatial model we must derive the total intracellular-amount dynamics by integrating over the cell volume. At early times in low temperatures (no metabolism) we have the following system for drug concentration and transport into a single cell: To define the total uptake rate for the entire cell we integrate the intracellular concentration dynamics with respect to the volume of the cell, 3 Translating the multiscale model to a simple continuum model

Optimisation of effective parameters for the simple continuum model
The average steady state profiles for the full, multiscale models were acquired by extracting 8 1D radial profiles from cross-sections through the centre of the spheroid slice corresponding to the lines = 0, = 0, = and = − and calculating the mean. This method was initially validated by comparing the simple continuum model with the average radial profile of the full model with no intercellular space/zero porosity (i.e., the model used in Figure 4A), both using the default parameter set, i.e., " #•• = " and #•• = (see Figure S2A). In order to optimise the effective parameter values required to fit the simple continuum model to the cell-based models with inclusion of intercellular space, " #•• and #•• were varied by up to three orders of magnitude either side of the default dimensional value (e.g., " #•• = [1×10 -3 , 1×10 3 ]× " discretised over a log-scale for 51 points distributed within the interval) and the minimum residual error at steady state between the continuum and cell-based models was identified according to the following formula: where ˆ u¢ represents the continuum model output, ̅ˆt šš{ represents the average cell-based model output, is radial distance (discretised at every µm from 0 to 750 µm, = 1: 751) and * indicates the steady state. This process was repeated for both intercellular space geometries (wide and narrow) and a range of feasible permeability coefficients (corresponding to LogD7.4 = 1, 2, 3, 4, 5). For plotting purposes in Figure 5F, cell-based models were compared with the continuum model by making use of the standard R 2 error metric: (S78)