An analytical approach for quantifying the influence of nanoparticle polydispersity on cellular delivered dose

Nanoparticles provide a promising approach for the targeted delivery of therapeutic, diagnostic and imaging agents in the body. However, it is not yet fully understood how the physico-chemical properties of the nanoparticles influence cellular association and uptake. Cellular association experiments are routinely performed in an effort to determine how nanoparticle properties impact the rate of nanoparticle–cell association. To compare experiments in a meaningful manner, the association data must be normalized by the amount of nanoparticles that arrive at the cells, a measure referred to as the delivered dose. The delivered dose is calculated from a model of nanoparticle transport through fluid. A standard assumption is that all nanoparticles within the population are monodisperse, namely the nanoparticles have the same physico-chemical properties. We present a semi-analytic solution to a modified model of nanoparticle transport that allows for the nanoparticle population to be polydisperse. This solution allows us to efficiently analyse the influence of polydispersity on the delivered dose. Combining characterization data obtained from a range of commonly used nanoparticles and our model, we find that the delivered dose changes by more than a factor of 2 if realistic amounts of polydispersity are considered.

where k b is the Boltzmann constant, T is the temperature of the fluid, η is the dynamic viscosity of the fluid and d is the diameter of the particle. The velocity is given by where g is the gravitational acceleration constant, ρ is the density of the particle and ρ f is the density of the fluid. The boundary conditions corresponding to the experimental geometry of interest are D ∂P (x, t) ∂x x=L + V P (L, t) = 0, and D ∂P (x, t) ∂x x=0 + V P (0, t) = αV P (0, t), where α ≥ 0 represents the ability of the nanoparticles to associate with the cells. For a single nanoparticle, we assume a point source initial condition at x = x 0 , To solve Equation (1) subject to Equations (4)-(6) we first non-dimensionalise the system by introducing the variables χ = x L , Pe = V L D , and τ = tV L , where Pe is the Peclet number. The model therefore becomes P (χ, τ ) ∂τ = 1 Pe subject to the initial condition and the boundary conditions and We next introduce the new variable Equations (8)-(11) can be expressed in terms of this new variable, which gives rise to the governing equation with initial condition and the boundary conditions and To solve this, we follow a standard separation of variables approach and assume that functions X(χ) and T (τ ) exist such that Q(χ, τ ) = X(χ)T (τ ).
Equation (13) therefore implies that where λ 2 are the eigenvalues of the system. Note that the dash refers to an ordinary derivative. The solution for T (τ ) is straightforward and is The solution for X(χ) depends on the sign of the eigenvalues. Interestingly, depending on the value of α, the number of non-positive eigenvalues change. Assuming λ 2 = 0, Upon substitution into the boundary conditions, it can be seen that there are no non-trivial solutions that satisfy the boundary conditions unless in which case there is a single non-trivial solution and hence the leading eigenvalue is zero. If we instead assume λ 2 < 0, X(χ) = C 1 cosh(λχ) + C 2 sinh(λχ).
Again, there are no non-trivial solutions that satisfy the boundary conditions unless in which case there is one non-trivial solution, which gives rise to a negative leading eigenvalue that can be obtained from the solution to Finally, if λ 2 > 0, X(χ) = D 1 cos(λχ) + D 2 sin(λχ).
There are an infinite number of eigenvalues that satisfy the boundary conditions for this form of the solution, and the eigenvalues can be obtained from the solutions to There is exactly one root of Equation (26), λ n , on each interval [nπ, (n + 1)π], except for the case where where the leading eigenvalue is non-positive as discussed previously, and hence (26) has no solution on [0, π].
Hence there are three solution regimes for the model depending on the values of α and Pe: To obtain the values for the coefficients A n we use generalised Fourier Series and the initial condition. That is, .
P (x, t) describes the probability distribution of the position of a particle initially located at x = x 0 and we note that there is an implicit dependence on x 0 . However, we are typically interested in more than a single particle. Hence we consider that there is initially some distribution of particles with the same properties throughout the fluid. In the absence of additional information we make the assumption that the particles are distributed uniformly at random. The nanoparticle distribution throughout the fluid,P (x, t), is therefore given byP Furthermore, we are interested in calculating the delivered dose, U (t), which is the proportion of the initially administered dose that has arrived at the cell population. This can be calculated by evaluating the difference between the initial nanoparticle mass present in the fluid and the nanoparticle mass present in the fluid at time t. Specifically, We can calculate these integrals analytically.
where, for n ≥ 1, B n is defined in Equation (40) and where, for n ≥ 1, B n is defined in Equation (40) and Equation (38) describes the delivered dose for a nanoparticle population with a single set of physicochemical characteristics, that is, the same density and diameter. We are interested in the delivered dose across a polydisperse population of nanoparticles. If the nanoparticle diameters follow a size distribution S(η|θ) where θ are the shape parameters that govern the distribution then the delivered dose for this population, U p (t) is given by noting that the value of η influences V , D and Pe depending on the type of polydispersity considered. Note that this delivered dose corresponds to the number delivered dose.
For the mass delivered dose we add an additional weighting function that corresponds to the relative volume of each nanoparticle, and normalise by the total volume of the population, where d(η) is the diameter of the nanoparticle.

Method of solution
As the integral in Equation (46) cannot be solved for an arbitrary size distribution, we use a numerical approximation of the integral consisting of Equation (38) evaluated at ∆ θ values of the size distribution.
We truncate the infinite series in each solution after n max terms and note that the solution is insensitive to the inclusion of additional terms. To obtain the eigenvalues defined by the solutions to Equations (24) and (26) we use Matlab's fzero function, with the interval of possible function roots defined as previously. We note that in specific cases, corresponding to nanoparticles where sedimentation is significantly larger than diffusion, oscillations in the solution can appear at short times, due to the presence of Gibbs phenomena. In these cases, a numerical solution should be considered.

Experimental details
To obtain the scanning electron microscopy data we incubated 400 nm PEG (40kDa) @ mesoporous silica nanoparticles in Milli-Q in a 12 well plate containing a silica wafer. Particles composed of PEG were synthesised using a mesoporous silica template according to a previous published method [1,2], and the silica template was not removed from the nanoparticle. The nanoparticles were taken from a stock solution at 4.4×10 6 nanoparticles/µL, and 45 µL of this solution was added to 1455 µL of Milli-Q. Prior to incubation in the well, the nanoparticle/Milli-Q solution was vortexed and briefly sonicated to assure uniform distribution throughout the solution. Particles were undisturbed for 4 h and allowed to settle. After 4 h of incubation, fluid was removed from the well. The remaining silica wafer was then gold-sputtered and imaged using scanning electron microscopy (Philips XL30 FESEM) with a 15 kV beam. Images of the nanoparticles atop the silica wafer were captured at 12,000x magnification.
To obtain the size distribution of the nanoparticles we chose three representative SEM images, shown in Figures S1(a)-(c). The size of the nanoparticles was detected and measured automatically using Matlab's imfindcircles function with a minium radius of 7 pixels, a maximum radius of 30 pixels and a edge threshold detection value of 0.3. We overlay the detected circles in Figures S1(d)-(f). Red circles correspond to detected circles that are manually removed due to inaccuracy, while green circles correspond to detected circles that contribute to the measured size distribution (Main manuscript, Figure 5(b)). Finally, we fit a lognormal distribution to the histogram of measured nanoparticle diameters.  Figure 1: Scanning electron microscopy images of 400nm PEG @ mesoporous silica nanoparticles. Green circles represent nanoparticles detected and measured by our image analysis algorithm. Red circles correspond to detected nanoparticles that are manually removed, and hence not measured, as inspection revealed these to be false positives.