Revisiting cell–particle association in vitro: A quantitative method to compare particle performance

Nanoengineering has the potential to revolutionize medicine by designing drug delivery systems that are both efficacious and highly selective. Determination of the affinity between cell lines and nanoparticles is thus of central importance, both to enable comparison of particles and to facilitate prediction of in vivo response. Attempts to compare particle performance can be dominated by experimental artifacts (including settling effects) or variability in experimental protocol. Instead, qualitative methods are generally used, limiting the reusability of many studies. Herein, we introduce a mathematical model-based approach to quantify the affinity between a cell–particle pairing, independent of the aforementioned confounding artifacts. The analysis presented can serve as a quantitative metric of the stealth, fouling, and targeting performance of nanoengineered particles in vitro. We validate this approach using a newly created in vitro dataset, consisting of seven different disulfide-stabilized poly(methacrylic acid) particles ranging from ~100 to 1000 nm in diameter that were incubated with three different cell lines (HeLa, THP-1, and RAW 264.7). We further expanded this dataset through the inclusion of previously published data and use it to determine which of five mathematical models best describe cell–particle association. We subsequently use this model to perform a quantitative comparison of cell–particle association for cell–particle pairings in our dataset. This analysis reveals a more complex cell–particle association relationship than a simplistic interpretation of the data, which erroneously assigns high affinity for all cell lines examined to large particles. Finally, we provide an online tool (http://bionano.xyz/ estimator), which allows other researchers to easily apply this modeling approach to their experimental results.

with a newly formulated particle, and then to measure the cell-particle "association" -a term that encompasses both membrane-bound particles and those that have been internalized [1].
Unfortunately, at present it is essentially not possible to meaningfully compare results from distinct association incubation experiments, for three reasons. First, results from incubation experiments are generally only reported qualitatively, for instance, by reporting the median fluorescence intensity (MFI) of cells after incubation with fluorescently labeled particles. Detected fluorescence is highly dependent on the labeling strategy and particle [2], instrumentation settings, environment [3], and cell autofluorescence [4], all of which can substantially vary from experiment to experiment. Second, the amount of material that reaches the surface of cells can be significantly different from the amount of material administered, due to particle movement in solution (e.g., sedimentation and diffusion) prior to contact with the cell surface. The effect of this "particle dosimetry" can vary by six orders of magnitude for different particle systems [5] and is a confounding variable when performing in vitro adherent cell culture experiments. While the impact of dosimetry effects was initially questioned [6], there now exists strong evidence, both experimental [7][8][9][10] and theoretical [5,11,12], that these effects are significant. Dosimetry effects are especially of concern when seeking to compare particles of different sizes or classes of material (e.g., inorganic vs organic) [13]; however, they are underreported in the bio-nano literature and are often unaccounted for in experiments [14]. Third, there are significant variations in the protocol used for incubation experiments, including the incubation time (with 1-24 h being common), administered particle concentration (which can vary from single particles to millions per cell), method of particle administration [15], number of cells added, and dimensions and type of cell culture (e.g., static vs adherent culture). Even if all details of an experimental protocol are provided, there is no currently established methodology to combine this information in order to quantitatively compare disparate experiments. Here, we present an approach to solve these issues.
One potential solution to these issues is to interpret experimental data using a mathematical model that accounts for both dosimetry effects and the kinetics of cell-particle association. With an appropriate model, kinetic parameters (e.g., rate constants) can be determined by fitting the model to experimental data. These kinetic parameters can then serve as a protocol-independent, quantitative metric of cell-particle interaction (Fig. 1). This approach has seen wide success in other areas of science. For instance, the Hill-Langmuir equation is extensively used in biochemistry and pharmacology [16] to quantify the interaction strength between two molecules, and compartmental models are a staple of pharmacokinetics [17]. In either of these contexts, it would be unusual to publish raw data without also fitting a field-appropriate kinetic model to data and reporting kinetic parameters. For the field of bio-nano research, there is currently no such routinely used approach.
The adoption of analogous, quantitative approaches would help advance the field [18].
To pursue a kinetic modeling approach, a model of cell-particle association must be chosen. Previous work on modeling bio-nano interactions has implicitly assumed various different models of association, for instance cells as associating with all particles in contact with the cell surface [12], association proportional to concentration in contact with the cell surface [11], as a Langmuir isotherm model at equilibrium [19], or probabilistically associating [20]. However, an investigation into which kinetic model of association best fits experimental data has thus far been lacking. Crucially, any kinetic model of association must also account for in vitro dosimetry effects, or it will be unlikely to produce useful insight, as fitting data to the model will be confounded by dosimetry effects.
Here, we compare five candidate kinetic models of cell-particle association, including some that have been previously described in literature, to determine what kinetic model of cell-particle association is most appropriate. To this effect, we synthesize a library of disulfidestabilized poly(methacrylic acid) (PMA SH ) particles ranging from~100 to 1000 nm in diameter and measure their association over time with three cell lines (HeLa, THP-1, and RAW 264.7 (subsequently referred to as RAW)) in separate time course incubation experiments. We further expand this dataset with previously published association experiments. We evaluate the performance of each kinetic model using all experiments from this dataset and then perform further analysis using the model with closest alignment to experimental data. Additionally, we demonstrate the problems that arise when a simplistic approach is used to evaluate particle-cell interaction that does not account for dosimetry or association kinetics. Finally, we discuss the implications of these results for particle design and performance evaluation, as well as limitations and future directions.

Particle library synthesis and characterization
Disulfide-stabilized poly(methacrylic acid) (PMA SH ) core-shell particles (core-shells) were synthesized according to a previously described method [21]. Two different types of monodisperse nanoparticles were used as templates for the preparation of multilayered films across the nano-to micro-range: SiO 2 nanoparticles (NPs) with average diameters of 235 and 519 nm were acquired from a commercial source, whereas gold nanoparticles (AuNPs) with average diameters of 60 and 110 nm were prepared using an adaptation of the seed-mediated growth procedure described by Perrault and Chan [22]. AuNPs were determined to be monodisperse by UV-Vis spectroscopy and TEM (Fig.  S1, Fig. 2)-the AuNPs had relatively narrow size distribution profiles and high circularity (Fig. S2). This was further confirmed through Fig. 1. Summary of cell-particle kinetic modeling. Characterization information about cells and particles along with details of experimental protocol are used as parameters in a given mathematical model, which accounts for both dosimetry and the kinetics of cell-particle association. Ultimately, the fit kinetic parameters, which are the output of this modeling, can enable unbiased, quantitative comparison between in vitro experiments that vary in cell line, experimental protocol, and/ or particle type. dynamic light scattering (DLS) measurements (Fig. S3). Four polymer bilayers were deposited onto either SiO 2 or AuNPs, using hydrogen bonding between the protonated carboxyl groups of PMA SH and poly(Nvinyl pyrrolidone) (PVPON). Disulfide cross-links were formed after depositing the four PMA SH /PVPON bilayers to stabilize the PMA SH multilayer films. The sacrificial PVPON layers were subsequently released from the films by raising the pH, thereby yielding monodisperse PMA SH core-shells. The resultant PMA SH polymer films had thicknesses of between 6 and 12 nm, corresponding to 2-3 nm per layer (Table S1), which is consistent with the thickness of PMA layers in previously reported thin films [23].
This library of PMA SH particles was expanded through the inclusion (core-shell) or exclusion (capsule) of the core template particle. To obtain capsules, PMA SH /PVPON core-shells were exposed to hydrofluoric acid to etch the SiO 2 particles or potassium cyanide to etch the AuNP particles, followed by PVPON removal (Fig. S4). The diameter of the particles was determined by TEM and DLS measurements (Table  S2). The hydrodynamic radius of capsules was approximately two times larger than that of the respective core-shells, which we attribute to swelling in physiological pH from the deprotonation of the carboxylic groups of PMA [24]. This deprotonation also provides all of the core-shells and capsules with negative ξ-potentials (Table S2) [25]. In the remainder of this work, we refer to synthesized PMA SH particles (core-shells and capsules) as per their hydrodynamic diameter as obtained by DLS measurements, rounded to the nearest nm (i.e., 633, 282, 150, and 95 nm core-shells; and 1032, 480, and 214 nm capsules). We generated a dataset of in vitro cell-particle association data with which to compare candidate kinetic models of association. We first confirmed that particles were nontoxic to each cell line up to particleto-cell ratios of 10,000:1 (Supplementary S1, Fig. S5) [26]. Next, we incubated the synthesized PMA SH particles with HeLa (adherent), RAW 264.7 (adherent), and THP-1 (suspension) cells for 1, 2, 4, 8, 16, and 24 h at a particle-to-cell ratio of 100:1. After incubation, cells were washed and analyzed via flow cytometry. Cellular association was further qualitatively confirmed via confocal microscopy (Fig. S6). These flow cytometry experiments formed the basis of our time course incubation association dataset. Additionally, we included flow cytometry time course association data from two projects previously published by our laboratory. The first dataset [8] involved two PMA particle systems (core-shell or capsule) incubated with adherent NIH/3 T3 fibroblasts in either static or continuously mixed incubation conditions, thus generating a total of 4 incubation experiments. The second dataset [27] consisted of 1.3 μm "targeted" hyaluronic acid-poly(ethylene glycol) metal-phenolic network (HA-PEG-MPN) particles. These particles were incubated either with a target cell line (adherent MDA-MB-231), expressing CD44 (to which hyaluronic acid binds), or a non-target cell line (adherent BT-474), thus generating a total of 2 incubation experiments. In total, our incubation dataset consisted of 25 different combinations of cell line, particle, and incubation type.

Simplistic comparison of in vitro data incorrectly assigns high affinity between cells and large particles
We report the two most common metrics (percentage cell association and MFI) used to analyze cell-particle association data in Fig. 3. This figure contains two separate analyses of all in vitro flow cytometry data for our synthesized PMA SH particles. Percentage cell association ( Fig. 3a-f) estimates the percentage of cells that have associated with a particle. MFI (Fig. 3g-l) reports the fluorescence intensity of cells and particles that have been incubated together (reported in this figure as MFI over background). Both these metrics indicate that cell association is greater for the larger particles, irrespective of the cell type examined (Fig. 3m). If we assume that these metrics accurately reflect the underlying biology of cell-particle interaction, this is a surprising result. As a particle increases in size, it is excluded from certain uptake pathways, so a priori we would expect smaller particles to exhibit higher association to cells. We offer this simple explanation for these seemingly anomalous results: larger particles are brighter than smaller particles, and larger particles settle out of solution faster than smaller particles. Thus, we conclude that percentage cell association and MFI do not fully reflect the underlying biology when comparing different types of particles, despite their wide usage within the published bio-nano literature.

A mathematical framework for dosimetric-kinetic modeling
It is important to distinguish between the association model, which describes the biological interaction between cells and particles, and the dosimetric model, which describes the movement of particles in solution prior to contact with cells. The model of association represents a high-level view of the biology and it is from this model that useful kinetic parameters may be determined. In contrast, the dosimetry model represents an experimental artifact of in vitro culture experiments; these effects must be accounted for but are not part of the underlying biology. We introduce here a mathematical framework that allows these two distinct parts to be composed together.
As in previous work on cell-particle dosimetry [11,12,28], we model particle concentration as a continuous and conserved quantity in solution and simulate particle transport through a fluid using partial differential equations (PDEs). In the most common case of particle incubation with adherent cells on the bottom of a plate well, cells form one boundary, open air forms the top boundary, and the plastic walls of the well form the remaining boundaries. Boundary conditions of a PDE simulation are naturally expressed as flux out of solution, i.e. rate of change of concentration over unit area ( u n d d ). In contrast, chemical reaction kinetics are typically both expressed and measured in terms of change of concentration over time ( u t d d ). It is thus necessary to establish a relationship between these two representations. Inspired by work drawing parallels between the kinetics of endocytosis and the law of mass action [29], as well as a result from Erban and Chapman [30], we accomplish this by considering cell-particle association f cell (mol m −2 s −1 ) as a reaction that occurs only on the "cell boundary" dΩ cell . This is represented by the following equation: where f fluid (mol m −2 s −1 ) is a function which represents particle motion in fluid (i.e., dosimetry); δ b (m −1 ) is a Dirac delta-like function with an impulse along the cell boundary surface dΩ cell , and is zero away from the cell boundary [31]; and dΩ are all boundaries of the solution. It can be shown that this is equivalent to the equation (full derivation in Supplementary S2): This formulation allows for the mathematically rigorous composition of a particle fluid transport model f fluid with any kinetic model of cell-particle interaction f cell . Thus, we expect this result to be useful in the development of kinetic models of cell interaction.
Here, we chose f fluid according to the type of incubation experiment. For conventional in vitro incubation of particles with adherent cells on the bottom of a plate in a static configuration, we choose f fluid to account for sedimentation and diffusion: where D (m 2 s −1 ) is the diffusion coefficient for particles in media (calculated using the Stokes-Einstein equation), and s (m s −1 ) is the sedimentation/advection vector (calculated using Stokes' law). This choice reflects the significant body of research that demonstrates that diffusive and sedimentary forces are responsible for variations in particle transport prior to contact with cells [5,7,8], as well as the relatively monodisperse nature of the studied particles (Fig. S3). In addition, we modeled the initial condition of our particle-cell solutions as being well mixed, consistent with our experimental procedure. We also use data from incubation experiments with suspension cells and with adherent cells under continuous mixing. We model both conditions as remaining well mixed for the entirety of the experiment. This simplifies the model to the following ordinary differential equation (ODE) for u(t), the concentration of particles at the cell boundary (derived in Supplementary S3): where S (m 2 ) is the surface area of the cell boundary and V (m 3 ) is the volume of the solution. This formulation allows for any choice of particle transport model to be used while preserving the same biokinetic model of cell-particle association. For instance, for particle systems that experience significant aggregation within fluid, we would choose f fluid to include explicit representation of aggregation/agglomeration such as the In Vitro Sedimentation, Diffusion and Dosimetry (ISDD) [5] model or the ISD3 [12] model. The dosimetric model f fluid could include convection currents within in vitro experiments, more advanced diffusion models, or more exotic incubation conditions such as cells under a . Regardless of analysis method, type of particle, or cell line, the overall trend is (m) a spurious correlation between particle size and association with cells, which can be explained by larger particles being brighter and/or settling faster than smaller particles.
continuous supply of flowing particles.
Unlike dosimetric effects (f fluid ), which have been studied in greater depth, it is not clear a priori what mathematical function should be chosen to represent the kinetics of cell-particle association (f cell ). Therefore, here we will consider five different functions based on different assumptions and determine which best fits the data. Each function represents a different hypothesis about which factors are important in cell particle association. Each function can be mathematically derived from the most complex kinetic function (the "full" model, below) under specific simplifying assumptions. We describe each function below.

Perfectly absorbing (PA)
Here, all particles that reach the cell boundary instantly associate with cells, regardless of concentration (i.e., mathematically, the cell boundary is perfectly absorbing). Note that "perfectly absorbing" refers to the boundary condition, which is distinct from chemical absorption (i.e., adsorption). This model has been previously used to investigate and explain differences in cell association due to in vitro dosimetric variation between particle systems [5,13]. This model is equivalent to the full model under the assumption that cell carrying capacity and surface particle capacity are never reached (i.e., are arbitrarily high) and that association rate r is arbitrarily high.

Linear (LIN)
where SC (dimensionless) is the surface coverage of the cells, r (m s −1 ) is the rate of cell association (i.e. the affinity between a particular cell-particle pair), and u (mol m −3 ) is particle concentration. In this model, particles associate with cells at a rate that is directly proportional to their concentration at the cell boundary. This model is equivalent to the full model under the assumption that cell carrying capacity and surface particle capacity are never reached (i.e., are arbitrarily high).

Surface flux (SF)
where S capacity (mol m −3 ) is the capacity for particles of the surface of the cells. In this model, particles associate with cells at a rate proportional to their concentration at the cell surface (as in the linear model), but the concentration and rate of uptake saturate as the concentration at the cell surface increases. This model is analogous to a Langmuir isotherm or Hill equation model, with "rate of association" being the dependent variable.

Cell carrying capacity (CCC)
where SC (dimensionless) is the surface coverage of cells, P capacity (mol m −3 ) is the carrying capacity of cells for particles, and P assoc (mol m −3 ) is the current number of associated particles. In this model, cells associate with particles until they reach the association capacity P capacity . As P capacity is approached, the rate of association approaches zero.
This model represents a combination of the surface capacity and cell carrying capacity models.
In the remainder of the manuscript, for brevity, we refer to these models as PA (perfectly absorbing), LIN (linear), SF (surface flux), CCC (cell carrying capacity), or FULL.

Comparison of kinetic model performance in vitro
Next, we sought to determine which of the above functions best describes the kinetics of cell-particle association. To address this question, we fit each combined dosimetric-kinetic model to data from each of the experiments in our dataset. Simulating any of these combined models requires knowledge of a number of parameters, nearly all of which are determined directly from the physicochemical properties of the particle (size, density), the cells (cell surface area, number of seeded cells), or from details of the experimental protocol (plate well dimensions, amount of media added, media viscosity, media density, temperature). Full details on how these fixed parameters are determined and used is provided in 'Materials and methods, Parameter estimation'. Flow cytometry data was used to estimate the average number of particles associated with each cell at each timepoint ('Materials and methods, Particle counting'). All models except for PA include r, an unknown "affinity" or "rate" parameter with units of m s −1 representing the interaction strength between a specific particle-cell pairing. This rate parameter was determined by fitting each model individually to experimental data. The CCC and FULL models include P capacity , the maximum number of particles that associate with a cell. This was determined heuristically from the asymptotic behavior of the experimental data ('Materials and methods, Solving and fitting kinetic models').
The quality of model fit to the data was compared by relative model performance p model , using the following equation: where err best is the lowest mean squared error of all models within that experiment and err model is the mean squared error between the model and experimental data. Model performance p model thus ranges from 0 (worst) to 1 (best) for each experiment. Relative model performance of the 5 models for all 25 experiments is visualized in Fig. 4. PA (purple) had the worst performance among all models considered across all experiments. CCC (blue) and FULL (green) had the best performance across the surveyed experiments. Additionally, there is no appreciable difference between the performance of the CCC and FULL models or between the SF and LIN models. This indicates that surface capacity S capacity does not play a significant role in the quality of model fit to data and is not a relevant phenomenon when considering cell-particle association kinetics at the doses considered herein. As the CCC model is both simpler than the FULL model and fits data equally well, the CCC model was chosen for subsequent kinetic analysis.
Though Fig. 4 shows that the CCC kinetic model is the best model among those considered, it does not show how well this model fits experimental data. We therefore plot each fit of the CCC model against experimental data in Fig. 5. In general, these results qualitatively indicate good fit across the wide variety of considered particles, cell lines, and experimental conditions, suggesting that the processes represented in the CCC model are sufficient to explain particle association kinetics. Quantitatively, the mean squared error of each model to each experiment are provided at https://figshare.com/articles/Fit_rates_and_ errors/7666058.

Effect of size, cell line, and particle type on association kinetics
Determination of association rate parameter r for each set of time course data allows for direct, quantitative comparison across these different in vitro experiments. We calculate the association rate parameter r for the CCC model by fitting it to experimental data, and the relationship between diameter and association rate for different sets of experiments is displayed in Fig. 6. All the results in this figure use the same scale and can be directly compared to each other. We emphasize that without a kinetics-focused approach, which also accounts for dosimetry, any attempt to compare in vitro association experiments across cell lines, particle classes, or experimental protocols has questionable validity. For example, simplistic interpretation of the raw cell-particle association data would lead one to conclude that all cells have a higher association rate for both larger and denser particles (Fig. 3m). In contrast, when confounding dosimetric factors, variations in experimental protocol, and particle fluorescence intensity are controlled for, a different picture of cell particle interaction is revealed. Cell line rather than particle size had the most obvious effect on cell-particle association rate, with macrophage-like (RAW 264.7) and monocyte-like (THP-1) cell lines demonstrating the highest association rates, consistent with their role in foreign body recognition (Fig. 6a, b). The highest r observed is between RAW 264.7 cells and the smallest capsule particles (> 4.9 μm s −1 and > 3.0 μm s −1 for 214 nm and 480 nm capsules, respectively) in which we were only able to provide a lower bound (see 'Materials and methods, Solving and fitting kinetic models'). Interestingly, RAW 264.7 cells demonstrated sharp differences for similarly sized core-shell particles, with r values~4 orders of magnitude lower (1.5, 1.3, and 0.6 nm s −1 for 150, 282, and 633 nm core-shells, respectively). This higher r value between RAW 264.7 cells and capsules sharply decreased for the largest capsules surveyed (3.0 nm s −1 for 1032 nm capsules). One explanation for these observations is that smaller particles and more flexible capsules are able to access a wider range of uptake pathways than the larger or less-flexible particles. THP-1 cells displayed a slightly higher r for the core-shell particles, but overall have remarkably consistent r (within 1 order of magnitude) for all particles surveyed. HeLa cells demonstrated slightly enhanced r for capsules, which persisted across all sizes surveyed (Fig. 6c). Additionally, our dataset contained a targeted particle (HA-PEG-MPN) cultured with its target cell line (MDA-MB-231) and a non-target control cell line (BT-274) (Fig. 6d). The association rate r between the particle and the target cell line was 2-3 orders of magnitude higher than that between the particle and non-target cell line, demonstrating genuine targeting behavior. However, we note that this r is not significantly higher than r between the largest particles and the phagocytic cell lines surveyed (Fig. 6d vs Fig. 6a, b). This suggests that avoidance of recognition by non-target phagocytic cell lines remains an important consideration for engineered particles.
Finally, our dataset contained PMA core-shell particles and capsules made of the same material and of approximately equal size, deliberately engineered to be as similar as possible [8,13], that were incubated with NIH/3 T3 fibroblasts in one of two experimental conditions: continuous mixing, or standard, static adherent culture (Fig. 6e). Thus, though these particles experienced different dosimetric effects and raw association data, they were expected to experience similar cell association kinetics. Encouragingly, the calculated r value between this cell line and the synthesized particles remained within a single order of magnitude, providing further support to the idea that estimating kinetic parameters provides a measure of cell-particle interaction that is largely independent of experimental protocol.

Conclusions
Here, we have combined particle synthesis, in vitro experiments, and mathematical modeling to present a model-based approach to analyzing cell-particle interactions, focused on determining the kinetics of cell-particle association. Our mathematical framework allows for the composition of a kinetic model of association, while accounting for dosimetric effects, which have been shown to be a significant confounder of in vitro incubation experiments. We compared the performance of five different kinetic models of association and found that a "cell carrying capacity" model was the simplest model to give good performance across surveyed experiments.
Analogous to kinetic modeling in other domains, the end product of this approach is a metric, in this case a rate parameter r, which represents the interaction strength between a particular pairing of cell line and particle. This can be used as a quantitative metric of cell-particle affinity and allows the comparison of experiments conducted using different materials, cell lines, and conditions. We report the rate parameters estimated across a range of materials and cell lines and find that estimated rate parameters are relatively insensitive to alterations in experimental conditions (well-mixed vs static, concentration variation for nonspecific interactions), supporting the claim that they can act as an experiment-independent quantitative metric of cell-particle Kinetic modeling and quantitative measurements of material performance remain a rich area of future research. We envision that potential stealth or targeted materials will one day be challenged against wide panels of cell lines, representing not only the target tissue, but also highly phagocytic "interfering" cell lines including macrophages and monocytes. In conjunction with the modeling framework we present here, the full targeting or stealth performance of a particle could be assessed in vitro, making the impact of particle engineering choices quantitative and explicit. Developing additional kinetic models also represents a promising direction of research. For instance, the dosimetric component could be expanded to incorporate non-spherical particles or alternative fluidic experimental conditions. The association component could incorporate cell-cycle variations or known surface receptor concentrations. A dataset that distinguished membrane bound and internalized particles or other sub-cellular localization data could lead to further development of multi-compartment kinetic models [28]. Additionally, as our focus has been on understanding the fundamentals of bio-nano interactions, within this work we have focused on association rate r as a function of particle number. However, other metrics of concentration [32] may be more appropriate for different research domains -for instance, mass, volume, or surface area concentration metrics may be more appropriate for more pharmacologically focused work.
We emphasize that the technique we present is inferential, rather than predictive: it enables quantification of in vitro experiments that were previously only qualitative. This quantitative approach enables the comparison of cellular response when particle physicochemical properties are varied (including surface modifications, charge, size, or stiffness). Modeling of this form can be easily incorporated as part of a high-throughput screening workflow, to enable rapid assessment of engineered particle systems. It is our hope that the tools and techniques we present will support quantitative and systematic comparison, reporting, and development of nanoengineered materials for biological applications.
Taken together, the present work represents a methodological "road map" for quantitative in vitro measurement and reporting of cell-particle association experiments. To facilitate this process for other researchers, we have developed a companion website located at http:// bionano.xyz/estimator, through which our modeling analysis can be applied to cell-particle time course association data following the MIRIBEL reporting standard [33]. Images of this site are displayed in Fig. 7. Furthermore, for research groups that are interested in developing other kinetic or dosimetric models, we provide a mathematical framework that allows for the composition and numerical solution of any kinetic model of association with any model of dosimetry, and we include a fully annotated version of our dataset to enable further analysis and development of kinetic models.

Particle synthesis
Prior to use, poly(methacrylic acid) conjugated to pyridine dithiotehylamine (PMA PDA ) was dissolved at a concentration of 100 mg mL −1 with 0.5 M dithiothreitol (DTT) solution in 3-(N-morpholino)propanesulfonic acid (MOPS) buffer (20 mM, pH 8) for at least 15 min at 37°C to expose free thiol groups (PMA SH ). Subsequently, the polymer solution was diluted with NaOAc (50 mM, pH 4) buffer to a concentration of 1 mg mL −1 just before layering. PVPON was directly dissolved in NaOAC (50 mM, pH 4) to a concentration of 1 mg mL −1 . In a typical experiment, a suspension of SiO 2 particles (100 μL, 5 wt% suspension) was washed three times with NaOAc buffer (50 mM, pH 4). Then, 500 μL PVPON (1 mg mL −1 ) in NaOAc buffer was added and incubated for 15 min with constant shaking. The PVPON-coated SiO 2 particles were then washed three times with NaOAc buffer and incubated with 500 μL PMA SH (1 mg mL −1 ). The suspension was incubated for 15 min and washed three times in NaOAc buffer. The outlined procedure describes the assembly of a single polymer bilayer. The adsorption of subsequent interacting polymers (PVPON/PMA SH ) was repeated until four bilayers were deposited. A final capping layer of PVPON was layered, resulting in ((PVPON/PMA SH ) 4 /PVPON) core-shells. The multilayered polymer film was then stabilized using thioldisulfide chemistry. Disulfide-stabilized core-shells were obtained by treating the suspension with 28 mM chloramine-T hydrate (CaT) solution in 2-(N-morpholino)ethanesulfonic acid hydrate (MES) buffer (50 mM, pH 6) for 1 min, followed by two washing steps with MES buffer and redispersion in NaOAc (50 mM, pH 4).
For layering on AuNPs, a similar procedure was followed: 80 nM of AuNPs was washed three times with purified water. Then, 500 μL PMA SH (1 mg mL −1 ) in NaOAc buffer (5 mM, pH 4) was added and incubated for 15 min with constant shaking. The PMA SH -coated AuNPs  were then washed three times with NaOAc buffer and incubated with 500 μL PVPON (1 mg mL −1 ). The suspension was incubated for 15 min and washed three times in NaOAc buffer. The adsorption of subsequent interacting polymers was repeated until four bilayers were deposited, resulting in (PMA SH /PVPON) 4 core-shells. Disulfide-stabilized core-shells were obtained by treating the suspension with 28 mM CaT solution in MES buffer (10 mM, pH 6) for 1 min, followed by two washing steps with MES buffer and redispersion in NaOAc (5 mM, pH 4).
The capsules were then formed by dissolving the SiO 2 templates using NH 4 F (13.3 M)-buffered HF (5 M) at pH 5 (Caution! HF is highly toxic and great care must be taken during handling) or dissolving AuNP templates using KCN (10 mg mL −1 ). The resulting capsules were washed twice with purified water then three times in Dulbecco's phosphate buffered saline (DPBS).
Additional details of synthesis and labeling are provided in Supplementary S4

In vitro association studies
HeLa, THP-1, or RAW 264.7 cells were seeded in a 12-well plate (1 × 10 5 cells per well) in 1 mL growth media with the addition of 10% (v/v) fetal bovine serum (FBS) and incubated with Alexa Fluor 647 (AF647)-labeled PMA SH core-shells or capsules at a particle-to-cell ratio of 100:1 with incubation time varied from 1 to 24 h at 37°C in a 5% CO 2 humidified atmosphere. After incubation, unassociated particles were removed from adherent cells by gently washing with DPBS. Cells were removed from plates by treatment with 0.25% trypsin solution (200 μL/well) for 5 min at 37°C, 5% CO 2 . Complete growth media (300 μL/well) was then added to inhibit trypsin activity. Cell suspensions were collected and centrifuged at 400g for 5 min. Non-adherent cells were collected and washed with DPBS three times via centrifugation at 400g for 5 min. The resulting cell pellets were resuspended in DPBS. Cells associated with core-shells and capsules were then determined through flow cytometry (Apogee microflow cytometer) by the acquisition of the signal from AF647. At least 10,000 cells were analyzed for each sample.

Particle counting
To compare model to experiment and fit kinetic parameters from experimental data, it is necessary to quantify particles associated per cell (P assoc (t)) from both experimental data and from our mathematical model. Particles associated per cell were estimated from fluorescent flow data using the following expression: real bio pmt exp background particle assoc (11) where P real_assoc (t) is the number of particles associated per cell at timepoint t, M exp is the median fluorescent intensity (MFI) of cells incubated with particles for time t, M background is the MFI of cells without particles, M particle is the MFI of particles in solution, C bio is a correction factor for changes in fluorescence between particles in solution and those associated with cells, and C pmt is a correction factor to account for different photomultiplier (PMT) voltage settings. While C bio can be estimated by integrating microscopy techniques with flow cytometry [34], all particles analyzed here used highly photo-stable fluorescent labels, and C bio was set to 1 for all systems. C pmt was determined experimentally by systematic variation of PMT voltage for the same sample (Supplementary S5).

Parameter estimation
Accurate experimental reporting and characterization are vital precursors to the use of any computational model [35]. The characterization parameters required to simulate and fit our models, and how we determined them, are outlined below and fully listed in Table 1. The exact parameters used in each simulation are also specified in the supplementary INI files, which fully describe an experiment (https://figshare.com/articles/FCS_and_INI_files/7623671).
In static adherent culture, we use a dosimetric model (Eq. 3) that requires a diffusion coefficient D as well as a sedimentation velocity s. Diffusion coefficient was determined using the Stokes-Einstein equation: where k B is the Boltzmann constant, T is temperature, η is dynamic viscosity (of the media), and r p is particle radius. T was 310.15 K (37°C) for all experiments and η was 0.00101 kg m −1 s −1 . Radius r p was dependent on the particle in question, and we used particle size as determined via DLS, as they characterize particle size in a hydrodynamic state (Table S1).
Sedimentation velocity s was determined via Stokes' law: where g is gravitational constant, ρ p is density of particle, ρ m is density of media, r is radius of the particle, and η is viscosity of media. ρ m was assumed to be 1000 kg m −3 . ρ p was determined in the following ways: For synthesized core-shell particles, ρ p was taken to be the density of the underlying template. These densities were as follows: 235 nm silica and 519 nm silica, 1850 kg m −3 ; 60 nm gold, 7853 kg m −3 ; and 110 nm gold, 12,554 kg m −3 . For synthesized capsules, density was estimated by assuming a PMA film thickness of~7 nm (Table S1, TEM measurements), and using this to estimate the portion of the capsule that was made of PMA (density = 1220 kg m −3 ). The rest of the capsule was assumed to have a density equal to that of the media (i.e., 1000 kg m −3 ). This resulted in the following densities: 214 nm capsules, 1041 kg m −3 ; 480 nm capsules 1019 kg m −3 ; and 1032 nm capsules, 1008 kg m −3 . Previously reported HA-PEG-MPN particles density was estimated in a similar way at 1001 kg m −3 , and previously reported PMA particles were determined to have densities of 1650 and 1060 kg m −3 for core-shell particles and capsules, respectively. The surface area of the cell boundary, S, was determined via the following equation: where c a is the surface area of a single cell and c n is the total number of cells. These were details particular to each experiment. S capacity was set to be the total surface area of cells divided by the cross-sectional area of the particles (i.e., the total amount of particles that could fix around the cells), according to the following equation: P capacity , the capacity of a cell line for a particle system, was determined from experimental data via the following heuristic: where d last is the derivative of the last two experimental points (measured in particles/min), and d max the maximum derivative between any two points. This heuristic is motivated by the behavior of the CCC and FULL kinetic models. In them, as time t → ∞, P assoc (t) → P capacity , assuming the supply of particles from solution is not exhausted. In all experimental data, the supplied particle rate (particle-to-cell ratio of 100:1) was much higher than actual cell association. Thus, if an incubation experiment has proceeded for sufficient time, the highest value of P assoc (t) is roughly equal to P capacity . Conversely, if time t is close to 0, the behavior of the kinetic model is roughly linear and P capacity can be set arbitrarily high without significantly altering association rate r. Thus, this heuristic determines if P capacity has been approached from the shape of the associated association data and sets this parameter accordingly.

Solving and fitting kinetic models
As the only way for particles to be removed from solution is by association with cells, we calculate the number of particles associated with a cell from our model in the following way: PDE models were solved in 1D by progressively refining a mesh until the grid convergence index [36] (GCI) proposed by Roache [37] between progressive mesh refinements was < 0.01, indicating convergence. Mesh size h y was initially chosen the following equation: The ODE form of our model was solved numerically using the python scipy [38] library, which provides a binding to the Fortran ODE solving routine lsoda. Association rate r was fit to experimental data for both PDE and ODE models through the use of Brent's method [39]. Brent's method is a root-finding algorithm for functions of one variable, which requires the underlying function to be repeatedly solved, and bisects a solution. Our implementation using Brent's method repeatedly numerically solved the PDE/ODE (calculated as described above), and the mean squared error between P real_assoc and calculated P model_assoc was used as the loss function for Brent's method.
In some cases (e.g., if particle transport to the cell surface is the limiting factor in cell uptake), only a lower bound on the rate constant r can be determined. Limits on r were determined by first simulating the PA model on the given experimental parameters and then fitting the LIN model to this simulated data.

Software implementation
All code used for this analysis is available at https://bitbucket.org/ mwfcomp/cell_particle_kinetic_models. Analysis was performed on files generated with the Flow Cytometry Standard (FCS). Experimental and characterization details were specified separately in INI files that were interpreted by accompanying code. Raw FCS data and INI files are available at https://figshare.com/articles/FCS_and_INI_files/7623671. Code written is primarily in python and makes use of the FEniCS [40], numpy [41], scipy [38], matplotlib [42], and pandas [43] software libraries.