On predicting heterogeneity in nanoparticle dosage

Nanoparticles are increasingly employed as a vehicle for the targeted delivery of therapeutics to specific cell types. However, much remains to be discovered about the fundamental biology that dictates the interactions between nanoparticles and cells. Accordingly, few nanoparticle-based targeted therapeutics have succeeded in clinical trials. One element that hinders our understanding of nanoparticle-cell interactions is the presence of heterogeneity in nanoparticle dosage data obtained from standard experiments. It is difficult to distinguish between heterogeneity that arises from stochasticity in nanoparticle behaviour, and that which arises from heterogeneity in the cell population. Mathematical investigations have revealed that both sources of heterogeneity contribute meaningfully to the heterogeneity in nanoparticle dosage. However, these investigations have relied on simplified models of nanoparticle internalisation. Here we present a stochastic mathematical model of nanoparticle internalisation that incorporates hitherto neglected biological phenomena such as multistage internalisation, cell division, asymmetric nanoparticle inheritance and nanoparticle saturation. Critically, our model provides information about nanoparticle dosage at an individual cell level. We perform model simulations to examine the influence of specific biological phenomena on the heterogeneity in nanoparticle dosage. Under certain modelling assumptions, we derive analytic approximations of the nanoparticle dosage distribution. We demonstrate that the analytic approximations are accurate, and show that nanoparticle dosage can be described by a Poisson mixture distribution with rate parameters that are a function of Beta-distributed random variables. We discuss the implications of the analytic results with respect to parameter estimation and model identifiability from standard experimental data. Finally, we highlight extensions and directions for future research. 2010 MSC 92B05


Introduction
either strongly bound to the cell surface or have been internalised by the cell are recorded [15,19].

Cell population model 124
To describe the behaviour of the adherent cell population, we consider a lattice-based random walk 125 on a two-dimensional square lattice of size X∆ by Y ∆ where ∆ is the lattice width [28] (Figure 1(d)).

126
The random walk is chosen to be an exclusion process, where each lattice site may be occupied by, 127 at most, one agent [28,29]. Here agents represent cells and ∆ is chosen to be equivalent to a cell Agents on the lattice undergo motility and proliferation events [31]. We make the assumption that 135 cell death does not play a significant role in the experiment, that is, that the nanoparticle dose is not 136 cytotoxic. For a motility event, an agent is randomly selected, with replacement, with probability M 137 during a timestep of duration τ . The agent, located at (x∆, y∆), attempts to move to a target site that 138 is randomly selected from the four nearest-neighbour sites located at (x∆ ± ∆, y∆) and (x∆, y∆ ± ∆).

139
The motility event will be successful if the target site is unoccupied, and will otherwise be aborted from an exponential distribution to a hypoexponential distribution, which has been shown to more 149 accurately reflect experimental observations [32]. Note that the number of stages in the proliferation 150 model do not necessarily correspond to the four stages in the cell cycle. In our model, an agent at 151 proliferation stage p is randomly selected, with replacement, with probability P p during a timestep of 152 duration τ . If p = Π, the agent simply progresses to the p + 1th stage of the cell cycle [32]. If p = Π 153 the agent, located at (x∆, y∆), attempts to place a daughter agent at a target site that is randomly 154 selected from the four nearest-neighbour sites located at (x∆ ± ∆, y∆) and (x∆, y∆ ± ∆). Similar 155 to motility events, a proliferation event will be successful if the target site is unoccupied, and will 156 otherwise be aborted (Figure 1(f)) [31]. If the event is successful, the proliferation stages of both the 157 parent and daughter agent are set to p = 1. If the event is aborted, the agent remains in stage Π. initially, on average, N initial /XY nanoparticles in the media above each lattice site (Figure 1(c)).

170
During a timestep of duration τ , an individual nanoparticle in the media can bind to an agent with probability P bind . This binding probability represents two processes. First, the nanoparticle must reach the cell layer during that timestep via random motion and an agent must be present on the relevant lattice site. Nanoparticle diffusivity arises from the Stokes-Einstein equation [26], and is given by where k b is the Boltzmann constant, T is the temperature of the media, η is the dynamic viscosity of the media and d is the nanoparticle diameter. The inverse relationship between nanoparticle diffusivity and diameter implies that smaller nanoparticles are more likely to arrive at the cell layer during a timestep. The second process is the nanoparticle binding to the cell, given that the nanoparticle has arrived at the cell layer. This occurs at a rate proportional to the strength of interaction, also known as the affinity, between the nanoparticle and the cell [19]. The probability that the nanoparticle arrives at the cell layer and binds to a cell during a timestep is a function of media depth, h, and nanoparticle diffusivity [35], P bind = λ 1 τ = β Dπτ 4h 2 , where 0 ≤ β ≤ 1 is the probability that the nanoparticle binds to the cell, given it has arrived at the cell 172 layer, and λ 1 is the rate parameter describing the transition from a nanoparticle in the fluid (s = 0) to 173 a nanoparticle bound to a cell membrane (s = 1). Calculating this probability allows us to determine 174 the rate of nanoparticle binding without explicitly tracking the locations of nanoparticles in the media.

176
The first stage of the internalisation process corresponds to nanoparticles binding to the cell mem-177 brane. Depending on the physicochemical properties of the nanoparticles and the cell phenotype, the 178 nanoparticle can be internalised via a number of endocytic pathways. These pathways can involve 179 macropinocytosis, receptor recruitment, formation of caveolae or membrane invaginations, and endo-180 somal trafficking [4,17]. Additionally, the pathway is not necessarily unidirectional, as nanoparticle-181 receptor binding is reversible and nanoparticle recycling (exocytosis) can occur. As, in general, these 182 pathways are not fully characterised, we choose to represent internalisation as an S-stage stochastic 183 process. We make the assumption that each stage of internalisation is a Poisson process, that is, 184 progressing through the stage requires an exponentially-distributed amount of time. We denote the 185 rate parameter of the sth stage as λ s , which represents the rate of transition from state s − 1 to state s. We perform all simulations of the model on a lattice with X = 10 and Y = 10 sites in the x and y directions, respectively. Initially, C 0 unique sites are selected at random and an agent is placed on each site. The cell cycle stage of each agent is selected uniformly at random. We use a fixed timestep approximation to simulate the evolution of the system in time [37]. In each timestep of duration τ we update the system according to three types of events: nanoparticle internalisation events, cell motility events, and cell proliferation events. For internalisation events, we keep track of the number of nanoparticles in each stage of internalisation for each agent in the simulation. We assume that nanoparticles are uniformly distributed throughout the culture media and hence, on average, progression to the next stage occurs with probability λ e s (jτ )τ . We select τ such that λ s τ 1 (and 213 M ≤ 1, P p ≤ 1 ∀ p), so that it is unlikely that a nanoparticle will progress through more than one 214 internalisation stage in a single timestep. We calculate the number of nanoparticles that progress 215 between internalisation stages for all agents and stages, and then simultaneously update the state of

230
We are interested in the distribution of N s (t) across the agent population, as this represents the 231 dosage distribution for each stage of the internalisation process. Typically, it is not possible to obtain 232 such information experimentally, as it is not straightforward to identify which stage of internalisation 233 a nanoparticle is in [15]. However, as we will demonstrate, it is possible to obtain analytic results for 234 the dosage distribution for each stage of internalisation under certain simplifying assumptions. The simplest form of the model that we consider is a single stage internalisation process for a cell population that is fully confluent. Here we will use the terminology of associated nanoparticles to highlight that there is no distinction between nanoparticles bound to the cell membrane and nanoparticles that have been internalised by the cell [19]. While this is an oversimplification of the internalisation process, it aligns with standard experimental practice [15]. As the cell population is fully confluent, no division events are possible, and nanoparticle inheritance can be neglected. We further assume that the impact of the carrying capacity can be neglected, that is, that K 1 → ∞ and hence λ e 1 (t) = λ 1 . Under these assumptions, the mean number of associated nanoparticles per cell by time t is which is the cumulative distribution function (CDF) of an exponential distribution with rate parameter λ 1 , multiplied by the average number of nanoparticles initially present in the culture media above 239 each lattice site. As this is simply an exponential process, the distribution of the number of associated  In mathematical models, it is standard to assume that experiments are performed on fully confluent cell populations; however, in certain cases it is desirable to commence experiments before the cell population becomes fully confluent [23]. The cell population may or may not reach full confluence by the end of the experiment, and this information is not commonly reported. If the cell population is not fully confluent, then the total number of associated nanoparticles is necessarily reduced as it is possible for nanoparticles to arrive at the base of the culture dish without interacting with a cell. If we assume that the amount of cell proliferation over the course of the experiment is negligible, then the number of associated particles per cell at time t is Poisson distributed with rate parameter where 0 ≤ C 0 ≤ XY is the number of cells in the population. This can be considered as the CDF of an If the cell population is not fully confluent initially, and cell proliferation cannot be neglected, it is not straightforward to obtain a closed-form analytic solution for the mean number of associated nanoparticles per cell. For a single stage model of cell proliferation (i.e. Π = 1), the mean number of cells in the population, C(t), is well-approximated by logistic growth, provided that the rate of cell motility is significantly higher than the rate of proliferation [31], and hence where ξ = P 1 /τ is the rate of cell proliferation. To the best of our knowledge, there is no analytic expression for the distribution of the number of cells in the population at time t. We can, however, use Equation (3) to find an integral expression for the mean number of associated nanoparticles per cell in the presence of cell proliferation Similar to Equation (2), this can be considered as the CDF of a random variable that represents the  distribution of associated nanoparticles per cell is noticeably broader if the inheritance is asymmetric (Figures 4(g)-(l) and 4(s)-(x)). In the asymmetric case, we observe that there is a substantial increase 292 in the number of cells that have zero nanoparticles associated by 24 hours, which only occurs rarely 293 for symmetric inheritance. We again observe that the dosage distribution is narrower for the multi-294 stage proliferation model (Figures 4(m)-(x)) than the single stage proliferation model (Figures 4(a)-(l)).

296
One benefit of employing the model is that it is possible to divide the population dosage distribution into separate dosage distributions based on the number of divisions a cell has experienced. That is, the population dosage distribution can be considered as a mixture distribution of the dosage distributions for cell subpopulations that have never divided, divided once, divided twice, and so on.
Here we make the assumption that the cell population is not close to confluence, and hence the growth of the cell population is approximately exponential. Further, we make the assumption that dosage depletion effects can be neglected, that is, that we are in the regime where association is independent of confluence. We first consider the subpopulation of cells that have experienced division exactly once by time t, which includes both the parent and daughter cells. We define the time at which the cell undergoes division as t 1 . Under the assumption of exponential growth, we know that t 1 is uniformly distributed on the finite interval 0 ≤ t 1 ≤ t. This can be seen by defining two i.i.d. exponential random variables X 1 and X 2 , which represent the time of the first and second divisions of a cell, and calculating p(X 1 = t 1 |X 1 < t, X 1 + X 2 > t) for 0 ≤ t 1 ≤ t. Accordingly, the dosage distribution of this subpopulation is a mixture distribution with two equally-weighted components: the dosage distribution of parent cells, which are Poisson distributed with rate parameter and the dosage distribution of the daughter cells, which are Poisson distributed with rate parameter recalling that I denotes the proportion of nanoparticles that are inherited by a daughter cell upon 297 division. There are two components for each rate parameter: the number of particles that associated 298 to the cell after division occurs, and the number of particles that associated before cell division, and 299 remained after division. There will be, on average, an equal number of daughter and parent cells that 300 have divided exactly once.

302
We follow a similar process to obtain the dosage distribution for the subpopulation of cells that have divided exactly twice by time t. We define the time at which the first and second divisions occur as t 1 and t 2 . From order statistics [38], we know that the PDF of t 1 is which implies that t 1 t ∼ Beta(1, 2).
In turn, the time between the first and second division events, t 2 − t 1 , is uniformly distributed on 0 ≤ t 2 − t 1 ≤ t − t 1 . There are four equally-weighted components in the subpopulation dosage mixture distribution, which are Poisson distributed with rate parameters We can generalise this process to obtain the dosage distribution for the subpopulation of cells that have divided N times by time t, where the time of the nth division (1 ≤ n ≤ N ) is denoted t n . Again, from order statistics [38], the PDF for t n − t n−1 (i.e. the time between division events) is noting that we define t 0 = 0, and hence To obtain the general form of the rate parameters, we define the sequence a n = {b 1 , b 2 , · · · , b n }, where b j ∈ {d,p}, which represents the order of the n division events that a cell has experienced, and whether it was the daughter (d) or parent (p) cell at each division. We define Ω d,i and Ω p,i as the number of times that a cell was the daughter cell or the parent cell, respectively, in the last i elements in the sequence. Accordingly, we define the rate parameter for an arbitrary sequence where t i for 1 ≤ i ≤ n is obtained from a Beta distribution as defined in (5) and t 0 = 0. We note that 311 there are 2 n permutations of the cell division sequence, and hence there are 2 n rate parameters for the 312 cell subpopulation that has experienced n divisions and the corresponding mixture distribution has 2 n 313 components.

315
We now highlight the validity of this approximate dosage distribution for each cell subpopulation.
All cells No division One division Two divisions Three divisions Four divisions  We next consider a two stage model, where nanoparticles first bind to a cell's membrane and, subsequently, become internalised by that cell. We assume that the cell population is fully confluent and that the carrying capacity at each stage of internalisation can be neglected. For a nanoparticle to be bound to a cell at time t, we require that the time of binding is less than t and that the time of internalisation is greater than t. If these two events occur according to two exponential random variables X 1 and X 2 with rates λ 1 and λ 2 , then the probability that an individual nanoparticle is bound at time t is equivalent to P (X 1 < t, X 1 + X 2 > t). Accordingly, the mean number of bound nanoparticles per cell is For a nanoparticle to be internalised at time t, we require that the combined binding and internalisation time is less than t, that is, P (X 1 + X 2 < t). Hence the mean number of internalised nanoparticles per cell is As for the single stage internalisation model, the effective rate parameter for binding decreases with where the cell population is fully confluent (Figure 7(a)) and where the cell population is not fully 373 confluent (Figure 7(b)). Further, these results highlight the importance of accurately distinguishing  It is possible to calculate the mean number of nanoparticles per cell for each stage of the internalisation process for an arbitrary number of stages in the absence of any carrying capacities. While we state the result here, we note that it would be exceedingly difficult to parameterise a model with more than two stages from the types of experimental data that are routinely collected, and hence we do not focus on analysing simulation results for S > 2 here. The number of nanoparticles that have reached the Sth stage of internalisation (i.e. that have been internalised) at time t is which arises from an S-stage hypoexponential distribution [39]. We can use this result to recursively define the number of nanoparticles in the sth stage of internalisation at time t for 1 ≤ s ≤ S − 1, as the number of nanoparticles in the sth stage is the difference between the number of nanoparticles that have reached the sth and (s + 1)th stages of internalisation We note that we can describe cell populations that are not fully confluent via a rescaling of λ 1 and the 382 number of nanoparticles available to each cell by XY /C 0 as previously. In the absence of a carrying If cell proliferation cannot be neglected, a Poisson dosage distribution is no longer appropriate. As in the case of the single stage internalisation model, it is possible to predict the dosage distribution by considering subpopulations of cells based on the number of division events that a cell has experienced. Again, we require that dosage depletion effects can be neglected, and that the initial confluency of the cell population is small. The distribution of cell division times is Beta distributed as defined in (5). For cells that have divided exactly once, the distribution of the number of bound nanoparticles per cell is a mixture distribution of equally-weighted two Poisson distributions with rate parameters is defined in Equation (7), and the superscript denotes that the rate parameters correspond to the first stage of the internalisation process. The distribution of the number of internalised nanoparticles per cell is a mixture distribution of two equally-weighted Poisson distributions with rate parameters where t 1 ∼ U (0, t), E[N 2 (t)] is defined in Equation (8), and the superscript denotes that the rate 387 parameters correspond to the second stage of the internalisation process.

389
We can express the components of the mixture distribution that represents the dosage of nanoparticles in the sth stage of internalisation for the subpopulation of cells that have experienced n divisions in general terms. Recall that a n = {b 1 , b 2 , · · · , b n }, where b j ∈ {d,p}, defines the number, order, and type of division events that a cell has experienced. The mixture distribution for the subpopulation of cells that have experienced exactly n divisions has 2 n equally-weighted components, each of which is a Poisson distribution with rate parameter where t i for 1 ≤ i ≤ n is obtained from a Beta distribution as defined in (5), t 0 = 0, E N s (t) is

423
Here we select K 1 = 10 and K 2 = 5, that is, at most ten nanoparticles can be bound to an individual  A key parameter that informs the efficacy of a potential nanoparticle-based therapeutic is the nanoparticle-cell affinity [19]. This parameter provides a quantitative estimate of the strength of interaction between nanoparticles and cells [19]. In the context of the model presented here, the nanoparticle-cell affinity is proportional to λ 1 . While this parameter can be estimated for each case where an analytic result is presented via numerical nonlinear curve fitting methods, provided that N s (t) is measured for all s, where the bar indicates that this is the mean of the experimental data, such methods can require specialist knowledge to implement. Although such an approach would be preferable, here we discuss heuristic approaches for estimating key biological parameters, which require less specialised knowledge. For example, from the result in Equation (1), we know that at early experimental time points we should observe that the rate of nanoparticle association is approximately linear, with the caveat that there is negligible cell division and a high carrying capacity. Accordingly, λ 1 can simply be estimated via the slope of nanoparticle association over time. This result can be extended for multistage internalisation processes, as λ 1 can still be estimated from nanoparticle association data. Provided that N s (t) is measured for all s and we have an estimate of λ 1 , other transition rates can be estimated heuristically from the data. For example, for a two-stage internalisation process, we can estimateλ

543
Appendix A. Experimental details 544 We aimed to assess the internalisation of 100 nm nanoparticles by Jurkat T cells. We hypothesised 545 that targeting the nanoparticles to different surface receptors would alter the uptake kinetics by trig-546 gering receptor-mediated endocytosis. To test this, we used a model nanoparticle system made of iron 547 oxide coated with protein G, which is straightforward to functionalise with a targeting antibody. The