## Abstract

In order to achieve selective targeting of affinity–ligand coated nanoparticles to the target tissue, it is essential to understand the key mechanisms that govern their capture by the target cell. Next-generation pharmacokinetic (PK) models that systematically account for proteomic and mechanical factors can accelerate the design, validation and translation of targeted nanocarriers (NCs) in the clinic. Towards this objective, we have developed a computational model to delineate the roles played by target protein expression and mechanical factors of the target cell membrane in determining the avidity of functionalized NCs to live cells. Model results show quantitative agreement with *in vivo* experiments when specific and non-specific contributions to NC binding are taken into account. The specific contributions are accounted for through extensive simulations of multivalent receptor–ligand interactions, membrane mechanics and entropic factors such as membrane undulations and receptor translation. The computed NC avidity is strongly dependent on ligand density, receptor expression, bending mechanics of the target cell membrane, as well as entropic factors associated with the membrane and the receptor motion. Our computational model can predict the *in vivo* targeting levels of the intracellular adhesion molecule-1 (ICAM1)-coated NCs targeted to the lung, heart, kidney, liver and spleen of mouse, when the contributions due to endothelial capture are accounted for. The effect of other cells (such as monocytes, etc.) do not improve the model predictions at steady state. We demonstrate the predictive utility of our model by predicting partitioning coefficients of functionalized NCs in mice and human tissues and report the statistical accuracy of our model predictions under different scenarios.

## 1. Introduction

Design and optimization of affinity–ligand functionalized nanocarriers (NCs) for diagnostic and therapeutic purposes remains an active area of biomedical research [1]. Immense advances have been made on the design front resulting in a huge repertoire of NC configurations that encompass carriers as diverse as liposomes, inorganic particles, DNA cages, nanodiamonds and polymer-based formulations such as polymerosomes and nanogels [2–4]. Each of these formulations shows higher efficacy only in a limited context and specific tissues, and in the majority of the cases the reason for the conducive behaviour is not always evident. This is due to a lack of understanding of the fundamental principles governing the pharmacokinetics (PK) of the NC [3,4]. The understanding of how the delivery vehicle interacts with the target tissue and how the design parameters may be varied to reach optimal function still remains preliminary, because of the large dimensions of the design space and the interconnected nature of the parameter space, as well as the complex physiological and biochemical factors that govern NC binding [3–5]. Traditional PK models, which are often guided by heuristic principles or empirically determined functions are inadequate in the context of targeted drug delivery, because their parameters can seldom be tuned rationally, and can only be obtained through exhaustive *in vivo* experiments. In contrast to the empirical approaches, biophysical approaches can enhance the predictive value of such models, especially for newer classes of vectors for which the design/parameter space is large/high-dimensional.

The need for rational design in biomedicine has long been recognized, and efforts to devise a framework to describe the complex behaviour of NCs have been taken both from the experimental and theoretical viewpoints. A combinatorial approach involves the synthesis of multiple NC formulations, each with well-controlled chemical and physical properties, and uses high-throughput screening techniques to map the relationship between the formulation-type and efficacy [6]. Theoretical studies take a systems-based approach which includes cooperative binding, the physical state of the cell membrane and entropy as synergistic/competing factors that determine NC avidity. Frenkel and co-workers [7,8] have used statistical mechanics-based models and Monte Carlo (MC) simulations to determine receptor–ligand interaction parameters for selective targeting of NCs. Weikl and co-workers [9–11] have also used a similar approach to investigate the effect of membrane roughness and lateral diffusion of receptors in mediating their interactions with ligands. Liu *et al.* [12,13] have developed a general computational framework that, in addition to explicitly representing the specific ligand–receptor interactions, also takes into account additional physiological factors such as the presence of glycocalyx barriers and shear forces due to blood flow. The internalization dynamics of NCs has also been explored in a number of recent studies [14,15].

Building on these theoretical developments that are faithful to the biophysical and thermodynamic factors, we propose to formulate a statistical mechanics-based model framework to predict the binding landscape of an NC in the native environment of a live endothelial cell in the vasculature, as well as to predict its tissue targeting behaviour across various organs in different species. We hypothesize that a comprehensive model for NC binding to live cells should account for a number of essential factors that include (i) the mechanical properties and protein expressions of the target cell membrane [10,16], (ii) surface functionalization of the NC [1,2,12,14,16], (iii) physiological and microenvironmental conditions such as blood flow, and margination due to the erythrocytes [2,16,17], (iv) specificity of receptor–ligand interaction [1,2,16,18], (v) non-specific interaction of the NCs in the target tissue [1,2] and (vi) receptor diffusion and translation and rotation of the NC [12,14]. In this sense, we present a next-generation framework that takes these various factors into account and predicts NC targeting to live cells in five different organs (lung, liver, heart, kidney and spleen) across two species (mouse and human).

The binding dynamics of an NC is strongly influenced by the specificity of the receptor–ligand interactions, molecular stiffness such as the flexure of proteins, deformations of the cell membrane expressing the target receptors, expression levels of the target proteins, physiologically mediated external forces, such as forces due to blood flow and red blood cell margination, and entropic forces due to thermal undulations. The multiscale computational platform presented in this article is shown in figure 1*a* and details are presented in §2. It is based upon the framework of equilibrium statistical mechanics and couples continuum field models for cell membranes with coarse-grained molecular scale models for the NC, antibodies and target receptors. The information flow between the molecular and macroscopic degrees of freedom is defined in terms of energy functionals (integrals) and interaction potentials, that are described in detail in §2. The conformational state of a functionalized NC interacting with the cell membrane is evolved through a set of seven types of Monte Carlo moves (see §2), which on sufficient sampling yields the equilibrium conformations of the combined system and the required statistics on multivalent binding. The model parameters illustrated in figure 1*a*, as inputs to the computational framework, are self-consistently determined either from experimental measurements or using molecular simulations which makes this framework to be a *zero-fit* model. The main output of the MC simulation is the calculation of the free energy landscape for carrier binding to cells, which is quantified through umbrella sampling and the weighted histogram analysis method. This free energy landscape is used to compute the avidity of NC binding.

The illustration shown in figure 1*a* highlights the major components of the proposed computational platform. These components can be broadly classified as: (i) a set of input parameters for the coarse-grained and continuum models that completely define the protein expression and mechanical properties of the target cell membrane, the biochemical interactions of the receptor–ligand bond and the flexural rigidity of the target receptors, and experimentally controllable quantities such as the geometry and the surface chemistry of the functionalized NC; (ii) a computational engine based on Metropolis Monte Carlo techniques to exactly compute the association constant *K*_{a}, for a specified mechano-chemical microenvironment and (iii) a framework that accounts for the targeted contributions due to NC binding to non-endothelial cells such as a macrophage. Most *in vivo* experiments only report the total concentration of NCs in the whole organ which also includes the binding and internalization of the carrier through a number of other non-specific pathways such as internalization through caveolin pits, and uptake by macrophages. Accurate representations of these non-specific contributions require a more detailed description of the tissue morphology and the physiological conditions. The modular design of our framework allows us to also include non-specific but highly relevant contributions such as physiological and microenvironmental conditions based on direct experimental measurements: these include non-specific capture mechanisms such as enhanced perfusion and oedema, capture through specific mechanisms not involving intracellular adhesion molecule-1 (ICAM1) (such as through the interaction of the antibody with FC receptors) [19], heterogeneity in tissue mechanics, and sub-cellular scale variations in the expression of the target proteins [2]. With these inputs from experiments, our model is able to predict the enhancement due to tissue targeting by quantitatively considering the biophysical factors outlined above. In the following sections, we describe the methodology and the procedure for parameter estimation before discussing the results.

## 2. Material and methods

### 2.1 Computational methods for nanocarrier binding to cell membrane

Monte Carlo protocols for the NC motion and adhesion to non-compliant surfaces have been extensively tested in previous works. Agrawal & Radhakrishnan [20] introduced the original version of the NC binding model using which they studied the effects of receptor flexure and the endothelial glycocalyx. Liu *et al.* [12] have extended this model by combining it with a powerful methodology for computing absolute binding free energies (described in the following subsection), through which they made successful comparisons to experiments probing NC binding to cells, tissue and NC unbinding in atomic force microscopy experiments. Later, Liu *et al.* [13] further extended their model to study the effects of shear flow and glycocalyx. Collectively, these works not only illustrate the model for NC binding in detail, but also explain how each parameter was obtained through independent well-controlled experiments, thereby obviating the need for any fitting of the data to obtain parameters. The values of some of the parameters determined from experiments (such as the flexural rigidity of the receptors) were also confirmed by independently carrying out molecular dynamics simulations [18]. The model presented in this work leverages the successes of the models described in the references cited above, and further extends these models by incorporating a crucial component, namely the compliance of the adhering interface.

The main components of the computational model are: (figure 1*b*) (i) the cell membrane, (ii) the position of the NC of radius *r*_{NC}, (iii) the coarse-grained positions of the *N*_{ab} antibodies defined on the surface of the NC and (iv) the coarse-grained positions and flexure of the *N*_{ant} receptors defined on the curvilinear manifold defined by the membrane surface.

#### 2.1.1 Model for the cell membrane

The patch of the cell membrane to which the NC binds is represented as a continuum surface whose curvilinear area is denoted by

For the purpose of numerical simulations, we represent the continuum membrane as a triangulated surface with *N*_{m} nodes, *T*_{m} triangles and *L*_{m} links, such that *N*_{m}+*T*_{m}−*L*_{m}=0. The discretized form of the membrane elastic energy [21,22] for the triangulated surface may be given in terms of the principal curvatures as [23]:

Here, *κ* and *σ* are, respectively, the effective bending rigidity and surface tension of the membrane of the target cell. *c*_{1,v} and *c*_{2,v} are the two principal curvatures at vertex *v* and *v*, such that the total curvilinear area is given by *a*_{0} and *a*_{0} is the discretization length-scale, which is taken to be *a*_{0}=10 nm.

#### 2.1.2 Model for the surface receptors

The receptor molecules are defined as cylindrical rods, with flexural angles *θ* and *ϕ* and length **a**_{b} and **a**_{t}, respectively. When *θ*=*ϕ*=0, i.e. in the unflexed state, the orientation of a receptor *i*, with flexural rigidity *κ*_{f}, is governed by the flexural energy:

#### 2.1.3 Model for the ligands

The ligand molecules are defined as rigid rods of length **A**_{b} are on the surface of the NC and their tip positions **A**_{t} are defined according to some predefined orientations—in the case of a spherical NC the ligand molecules are defined along the radial direction. The binding interaction between a receptor *i* and a ligand *j* is modelled using the Bell potential:
*d*_{ij}=|**A**_{t,i}−**a**_{t,j}| is the distance between the tip positions of the chosen receptor–ligand pair and *d** is the range of the binding interaction. *κ*_{b} is the spring constant of a receptor–ligand bond.

The total energy of the NC–membrane system is given by

#### 2.1.4 Monte Carlo moves

The conformational states of the system are evolved using a Metropolis Monte Carlo method that consists of seven independent moves: namely (1) a vertex move to simulate thermal fluctuations in the membrane, (2) a link flip to simulate membrane fluidity, (3) translation of the NC, (4) rotation of the NC, (5) diffusion of receptors, (6) receptor flexure move and (7) formation and breakage of receptor–ligand bonds. Moves (1–5) are performed according to the rules of canonical Monte Carlo while moves (6) and (7) are treated as rare events and are performed using the Rosenbluth sampling technique. Details of the various Monte Carlo moves are given in the electronic supplementary material, §S1. We note that with respect to the system (especially the membrane), we work in the ensemble of constant *T* is the system temperature.

### 2.2 Mean field model for cytoskeleton

The model for NC adhesion to the cell interface is expected to be influenced by the state of the cytoskeleton. While our model does not have an explicit representation for the cytoskeletal elements, the effects of the cytoskeleton are incorporated at the mean field level through the renormalization of the physical properties of the membrane, namely the values of *κ*, *σ* and *et al.* [25] have modelled a cytoskeletally fortified membrane using the Helfrich model by using a renormalized set of parameters described above for describing the mechanism of formation of the immunological synapse.

Several works in biomechanics have aimed to characterize cells based on mechanical measurements using a wide range of techniques such as flow and optical cytometry, manipulation using micropipette, optical tweezers and laser traps, and microfluidic devices (see [26–28] for comprehensive reviews). These studies have focused on whole-cell measurements. In the case of NC adhesion, the changes in mechanical properties are primarily caused by variations in the structure and organization of the cellular cytoskeleton [29] and the extracellular matrix [30]. Such sub-cellular scale rearrangements can significantly impact the mechanical properties of the cell membrane at length-scales smaller than cellular dimensions (i.e. tens of nanometres to less than 1 μm), a range which also corresponds to the scale at which the NC engages through multivalent interactions. The sub-cellular scale relevant to this discussion corresponds to the dimensions primarily set by the cortical cytoskeletal mesh, which has been estimated to be between *l*_{c}=150–500 nm [31,32]. The mechanical properties of a patch of the cell membrane that spans the region between multiple cytoskeletal pinning points, with typical dimensions *l*_{c}, can differ from the bulk because the nature of the thermal undulations (and the associated conformational entropy of the membrane) depends directly on *l*_{c}. In addition, the total area of the membrane

Drawing inspiration from the work of Qi *et al.* [25], we sought to describe the cytoskeletally fortified membrane using a renormalized set of parameters for the Helfrich model. We do so by introducing pinning sites on the membrane, which represent the binding of the cytoskeletal proteins to the membrane (figure 2). The pinning density is chosen to vary between 0 and 12% to mimic the distribution of *l*_{c} noted above [31,32]. We note that the lower end of the pinning density corresponds to *l*_{c}>500 nm and the upper end corresponds to *l*_{c}∼150 nm.

Even though the simulations are performed in curvilinear coordinates, for the sake of analysis of membrane undulations alone, we parametrize this surface in the Monge gauge and take the Cartesian *x*–*y* plane to be the reference surface—in this representation, every point on the membrane may be denoted by [*x*,*y*,*h*(*x*,*y*)] where *h*(*x*,*y*) is the height displacement along the *z*-direction. At non-zero temperatures, the spectrum of height undulations follow:
*q* is the wavenumber, *k*_{B} the Boltzmann constant and *T* the absolute temperature, where *h*_{q} is the Fourier transform of *h*(*x*,*y*).

The results for the height–height undulation spectra (see equation (2.5)) for membranes with pinned sites are compared with those in the absence of pinning (figure 2). The results show that for all excess areas considered, the effect of cytoskeletal pinning is to renormalize the values of *κ* and *σ* and that the scaling in the spectra can be well described by those predicted by the Helfrich model with renormalized parameters. In particular, the pinning can alter the value of *κ* in either direction depending on the *σ*. Our results show that for length-scales that are comparable to *l*_{c}, we can use the Helfrich model to investigate the binding free energy landscape for NC, and the renormalized parameters in figure 2 will enable us to consider effects of the cytoskeleton for a given state (represented by the pinning density and

### 2.3 Free energy analysis

#### 2.3.1 Umbrella sampling and weighted histogram analysis

The potential of mean force (PMF) for NC–membrane interactions is calculated using a macroscopic order parameter Δ*R*. This reaction coordinate is defined based on the combined conformational states of the NC and the membrane. If **R**_{CM} denotes the centre of mass of all membrane vertices within a distance of 2*r*_{NC} from the centre of the NC (denoted by **R**_{NC}) then the macroscopic order parameter is defined as the distance Δ*R*=|**R**_{CM}−**R**_{NC}|. Umbrella sampling using a biasing potential **U**_{bias}=*k*_{bias}(Δ*R*−Δ*R*_{0})^{2}/2 is performed at predefined values of Δ*R*_{0}, with a window interval of *δR*=2 nm. The strength of the biasing spring is chosen to be *k*_{bias}=2 *k*_{B}*T*/(*δR*)^{2}, i.e. the NC position can sample the entire window when it acquires 1 *k*_{B}*T* of thermal energy. Each biased window is independently sampled for 900 million Monte Carlo steps (one half of the moves are distributed equally between moves (1–5) and other half are distributed between moves (6) and (7) in the ratio 2 : 5) and the histogram of Δ*R* is recorded every 10 MC steps. The histograms of Δ*R* from multiple windows are combined and unbiased using the Weighted Histogram Analysis method in order to compute the PMF denoted by *R*, and hence is a true measure of the stability of such configurations. Here *β*=(*k*_{B}*T*)^{−1}. In our studies, we have used the thermodynamic integration (TI) technique to sample configuration windows that are prone to endpoint catastrophe and determine the additive constant to be the relative free energy computed using TI [33].

We run our code on 32-cores in parallel, and each PMF calculation (this includes TI and WHAM) takes 32×48 CPU hours. Each PMF is run in quadruplate (for computing standard deviations), and in total four values of *κ*, four values of *N*_{ant} and three *N*_{ab} densities were explored leading to an aggregate of 368 650 CPU hours of computing. The aggregate number of MC steps totalled 8640 trillion.

#### 2.3.2 Formulation of equations for binding avidity

In statistical mechanics, the binding avidity of a reaction is an exactly computable quantity and it has been applied with success in a number of problems involving small molecules and functionalized NCs [7,11,12,34,35]. Here, we present an equivalent form in terms of the PMF and the various entropic terms shown in figure 5 (details of the derivation are given in the electronic supplementary material, §S2). The association constant *K*_{a} for an NC functionalized with *N*_{ab}/NC antibodies, forming *n*_{b} simultaneous bonds upon binding to a membrane surface expressing *N*_{ant} receptors is given by

Here *L*_{z} is the dimension of the simulation box along the reaction coordinate, and *r** denotes the cutoff length at which the conformational states of the NC and those of the membrane–receptor system cease to overlap, i.e. the NC can only exist in an unbound state.

### 2.4 Pharmacokinetic model for tissue targeting

Experimental measurements of NC targeting in an organ are typically expressed in units of ‘percentage injected dose per gram of tissue’ (denoted as %idg) and are performed at the tissue scale, though the binding interactions for commonly used NCs (sizes in the range of 50–250 nm) occur at the scale of a single cell.

The association constant, given in equation (2.6), can be used within the framework of PK models to determine the effective partitioning coefficient of the targeted NC within the tissue at steady state. This modified framework (described in detail in the electronic supplementary material, §S3) explicitly takes into account (i) the non-specific binding of NC represented by *K*_{EC}, the association constant for the targeted adhesion of NCs to the endothelial cells and (iii) *K*_{M}, the association constant for an NC adhering to a macrophage. The standardized uptake value for the NC, with injected concentration *C*_{out}, in terms of *K*_{EC} and *K*_{M} is given by

Here, *ϕ*_{X} and *D*_{X} denote the concentration and diameter of endothelial cells (*X*=*EC*) and other cells (*X*=*M*) in the target tissue. The variable *r** defined in equation (2.6) for the different cell types and the variable *L*_{cap} represents the size of the cell free layer in the capillary in which the NC is perfused.

We present all %idg scores in scaled units given by *η*=(%*idg*)_{org,sp}/(%*idg*)_{lung,sp}, where the subscripts ‘org’ and ‘sp’ represent the target organ and species, respectively. The corresponding predictions from the computational model and experiments are denoted by *η*_{sim} and *η*_{exp}, respectively. We note that in these ratios the values of *L*_{cap} and *C*_{out} do not feature, while all other parameters are summarized in §2.7. We predict the endothelial targeting of anti-ICAM1-coated NCs in five different organs in mouse and compare them to targeting levels measured in *in vivo* experiments. The error bars in *η*_{sim} denote the standard deviation in the %idg scores, normalized by the %idg score for the lung. This standard deviation is determined by computing %idg scores for *K*_{EC}−Δ*K*_{EC}, *K*_{EC}, *K*_{EC}+Δ*K*_{EC}, where *K*_{EC} denotes the mean value of the association constant and Δ*K*_{EC} its standard deviation.

In our model, we note that the factor *In vivo* measurements (described in the next section) are performed in a timescale where the NC internalization will be minimal (as explained in that section). This justifies our neglect of NC internalization in the current formulation. We do include the fraction of NCs captured on the EC in the %idg calculations outlined above, which implies that NC fraction captured on the EC are accounted for in the tissue targeting estimates.

We further note that, for longer times, the effect of internalization can be included by augmenting the

### 2.5 *In vivo* targeting to vascular endothelium in mice

Anaesthetized C57BL/6 female mice (16–24 g, Harlan) were injected intravenously via jugular vein with NCs coated with murine anti-ICAM1 (YN1 clone, Biolegend) or control rat IgG (Jackson Labs). The injected dose was approximately 200 μl (or approx. 10 mg kg^{−1}) with a tracer amount of antibody-coated ^{125}I-labelled NC. Blood was collected from the retro-orbital plexus at 30 min post-injection and organs (heart, kidneys, liver, spleen and lungs) were collected at 30 min post-injection. Radioactivity and weight of the samples were determined to calculate NC targeting. These studies were carried out in accordance with the Guide for the Care and Use of Laboratory Animals as adopted and promulgated by the US National Institutes of Health. For a direct comparison, we quantified endothelial targeting as a function of organ uptake of NCs in mice. Full coverage of NCs was expressed as 100% endothelium targeting. We note that in the time of 30 min post-injection, the effect of ICAM1-mediated internalization of NC to tissue should be minimal because the timescale for such an internalization is about hour [36,37].

### 2.6 Comparison between model and *in vivo* experiment

We compare the %idg values computed from our model under various scenarios to those measured in experiments. In order to assess the predictive accuracy of our model, we compute the Pearson's correlation coefficients [38] (denoted by *r*^{2}) between model and experiment. The calculation of *r*^{2} considers both mean as well as the standard deviation in the computed and experimental values. We do so through a bootstrapping procedure, where five sets of 2000 datapoints are generated from a Gaussian distribution based on the mean and standard deviation of the calculated or experimental data. For each set of 2000 points, the *r*^{2} value is computed between model and experiment. We report the mean *r*^{2} over all 10 000 data points and a standard deviation in *r*^{2} based on the *r*^{2} values from each of the five sets.

The significance value (*p*-value) for the *r*^{2} for each model prediction is established by comparing the *r*^{2} for our tissue targeting model with that from a model representing the null hypothesis corresponding to no targeting. Specifically, the model corresponding to the null hypothesis assumes that the %idg is given by just *p*-value is computed using an unpaired *t*-test between our model and the model for the null hypothesis.

### 2.7 Parameter estimation

While our earlier works have discussed the methodologies for estimating model parameters [12,13,39], the detailed parameter sets used in our current model and the corresponding references are summarized for completeness in table 1.

#### 2.7.1 Estimating model parameters for ICAM1-targeted nanocarriers in mouse

In order to make reliable predictions for an NC binding to the target tissue, it is essential to accurately estimate the model parameters given in equation (2.7). This involves two set of measurements: (i) estimation of (*κ*, *N*_{ant})_{X} to compute *K*_{X} for *X*=*EC* and M, and (ii) estimates for *ϕ*_{X}, *D*_{X} to compute the biodistribution.

#### 2.7.2 Estimating protein expression levels

The expression levels of ICAM1 molecules (*N*_{ant}) for five major organs in mouse—namely the lung, liver, kidney, heart and spleen—are estimated by combining data from *in vitro* flow cytometry measurements [50], *in vivo* radio labelling measurements [50], mRNA expression levels reported in the BIOGPS database [51] and mass spectrometry measurements of proteome levels reported in the MOPED and PAXDb [52] databases. The number of ICAM1 molecules per square micrometre of the cell surface is shown in the main plot of figure 3*a*, and our results show that *N*_{ant} varies between 2 and 3 orders of magnitude across the various organs, with the expression being largest in the lungs (approx. 2000 ICAM1 μ*m*^{−2}) and smallest in the heart (approx. 50 ICAM1 μ*m*^{−2}), and the former is consistent with previous estimates for ICAM1 in mouse lungs [2]. The values of *N*_{ant} estimated from mRNA and whole-proteome measurements correlate very well with those measured using flow cytometry experiments, and this is shown in the inset to figure 3*a*. We note that while the correlation is good (with a reported *r*^{2}=0.98), the direct measures of accessible protein levels are the best choice for our model, where such data are available. Furthermore, we have also determined *N*_{ant} for other cells (e.g. resting monocytes and activated macrophages) using data from the BIOGPS database, and the corresponding values are given in figure 3*c*. Later in §3, we compute NC adhesion for the range of ICAM densities, *N*_{ant}=200, 500, 1000, 2000 and 4000 ICAM1 μ*m*^{−2}. We note that these densities are based on the projected area *N*_{ant} and

#### 2.7.3 Estimating mechanical properties of the target cell membrane

The bending rigidity and interfacial tension for cells are obtained from the work of Pontes *et al.* [57]. Based on our recent work [58], we have used a tether pulling assay to investigate the sub-cellular scale excess areas in different types of cells. Our estimates showed the bending rigidity and the excess area for most endothelial cells to be in the range of 40–60 *k*_{B}*T* and 10–30%, respectively. On the other hand, the bending rigidity of resting and activated macrophages varied substantially with *κ*=160 *k*_{B}*T* and *κ*=40 *k*_{B}*T*, respectively, while the excess areas for both these classes of macrophages were found to be around 20%.

The parameters for *κ* and *σ* obtained from the experiments noted above are already renormalized parameters including the effect of the cytoskeletal effects and can therefore be directly used in our Helfrich model. The length-scale of the membrane undulations probed in our model is set by the NC size which closely aligns with that explored experimentally by Pontes *et al.* [57] from which the representative values of the bending rigidity and membrane tension for the cell membrane patch are derived. Moreover, the

#### 2.7.4 Physiology-dependent parameters

The non-targeted partitioning coefficient *b* shows the values of *ϕ*_{EC}) and of macrophages (*ϕ*_{M}) in the target tissue varies across different organs, with the latter being high in clearance organs such as the liver and spleen and negligible in organs such as the lung and heart. Based on previously reported values [60], we set *ϕ*_{EC}=30% for all organs, and *ϕ*_{M}=3% for spleen and liver and 0% for the rest.

#### 2.7.5 Parameters for nanocarrier characterization

In table 2, we present the simulation parameters for NC characterization obtained from various sources in the literature.

#### 2.7.6 Interpolating the association constant *K*_{a}

An important question that arises in post-processing is how to estimate the association constant for arbitrary values of the receptor density, while we only have data for 200, 500, 1000, 2000 and 4000 ICAM1 μ*m*^{−2}.

A strategy we have adopted is to linearly interpolate the total free energy *β*=(*k*_{B}*T*)^{−1}. Consider two receptor densities *n*_{1} and *n*_{2} with corresponding free energies *K*_{a1} and *K*_{a2}. The free energy for an intermediate concentration *n** can be linearly interpolated to be
*K*_{a*} then takes the form:

## 3. Results and discussions

### 3.1 Can mechanical properties of the cell membrane be important determinants of nanocarrier avidity?

We hypothesize that (i) variations in the mechanical properties (*elastic parameters*) of the target cell membrane (the bending stiffness, *κ*, and the membrane excess area, *N*_{ant}) and (iii) the expression level of the affinity ligand on the NC (*N*_{ab}/NC) are the key determinants of NC avidity. We have tested these hypotheses by systematically varying *N*_{ab}/NC, *N*_{ant}, *κ* and

The equilibrium macroscopic conformations of a cell membrane are primarily governed by the two elastic parameters *κ* and *κ* varies between 20 and 200 *k*_{B}*T* [24] and the membrane excess area may vary from 0 to 300% [61,62], depending on the lipid composition and the degree of cytoskeletal pinning. In an earlier study, we have estimated this range to be between 1 and 35% [58] depending on the cell type. We note that the planar substrates used in earlier NC adhesion models characterize membranes with *κ* and *a* shows an anti-ICAM1-coated NC, with *N*_{ab}=162/NC, bound to cell membranes with varying values of *κ* and *κ*=20 *k*_{B}*T*) and stiff (*κ*=160 *k*_{B}*T*) membranes, with excess areas that fall into three distinct regimes: (1) small (*N*_{ant} (figure 4).

### 3.2 Multivalency of nanocarrier binding is a strong function of the mechanical variables of the binding surface

Multivalency, or the number of simultaneous (anti-ICAM1)–ICAM1 bonds formed between the ligands on the NC and the receptors expressed on the cell membrane is an important marker of NC avidity—it is a direct measure of the enthalpic contributions to the binding process. Figure 4*b* shows *P*(*n*_{b}), the probability distribution of the multivalency *n*_{b}, computed at equilibrium for an anti-ICAM1 functionalized NC in contact with each of the seven classes of cell membranes considered earlier; here *N*_{ab}=162/NC and *N*_{ant}=2000 ICAM1 μm^{−2}. The number of multivalent bonds formed by the NC is highly sensitive to variations in both *κ* and *P*(*n*_{b}), which for stiff membranes is normally distributed around a peak value of *n*_{b}≈8, for all the three cases of *n*_{b}≈10, for large values of *P*(*n*_{b}) for the membrane in the small deformation limit (panel (ii)) is consistent with that seen for flat substrates (panel (i)) and also agrees with that reported in earlier studies using planar substrates [12,13]. These results validate the hypothesis that *κ* and *n*_{b}) or if the entropic terms play a significant role as well. For example, it has been previously shown by us that in the case of flat substrates, the loss in the translational entropy of the diffusing receptors competes against the enthalpic gain due to receptor–ligand binding, and the interplay between these contributions limits the peak value of *n*_{b}≈3. In the following, we explore how such an effect manifests for cells with different membrane properties.

### 3.3 Bound receptors show different localization patterns under varying mechanical properties of the binding surface

The area traversed by a bound ICAM1 receptor on the cell surface is an indicator of the loss of translational entropy of the receptors on a cell surface due to multivalent engagement of receptor–ligand bonds [12]. Our results in figure 4*b* show the spatial map *P*(*t*_{1},*t*_{2}) to find a bound receptor–ligand pair at a point [*t*_{1},*t*_{2}] on a tangent plane defined by the mean orientation of the bound receptors (see the electronic supplementary material, §S2)—here the point [0,0] denotes the projection of the NC's centre of mass onto this tangent plane. It may be seen from figure 4*b*,*c* that: (1) receptors on a flat substrate that are bound to the NC are localized to an annulus-like region, consistent with that seen in earlier studies [12] and (2) in the presence of thermal undulations of the membrane (which is absent in the flat substrate model), *P*(*t*_{1},*t*_{2}) shows a larger spread and a diffuse pattern pointing to how bending modes in the membrane can manifest as additional binding modes for receptors on its surface. Our results, which display widely varying spatial maps of receptor–ligand pairs for different contexts of the cell membrane, clearly indicate that the translational entropy contribution can be cell-type (context) specific. Other enthalpic and entropic terms can also become significant in determining NC adhesion and these are discussed in the following.

### 3.4 Potential of mean force, which determines the enthalpic as well as entropic contributions to the binding avidity, is a strong function of the mechanical properties of the binding surface

The PMF *d*, and we find that both the form and the minimum value of the PMF (*κ* and

### 3.5 Expression levels of target proteins differentially affect nanocarrier targeting under varying mechanical properties of the binding surface

The density or the expression level of the target ICAM1 molecules is an obvious major factor that determines the binding affinity of the NC [2]. Modulations in ICAM1 density not only modulate the binding specificity of the NC, but also change its selectivity to the target cell [63]. We hypothesize that the effect of ICAM1 expression on NC binding will depend on the mechanical properties of the binding surface. In order to test this hypothesis, we systematically examined the various measures that quantify the binding characteristics of an anti-ICAM1-coated NC interacting with a cell membrane (with *κ*=40 *k*_{B}*T* and *N*_{ant}=200, 500, 1000, 2000 and 4000 ICAM1 μ*m*^{−2}, that are representative of the lung endothelium of a mouse with under-expressed, normal and over-expressed ICAM1 levels. In figure 5, as before, we quantify the various contributions to NC avidity using a set of measures, namely (1) *P*(*n*_{b}) (figure 5*a*), (2) *b*), (3) Δ*ϕ*Δ*θ*Δ*ψ*, the rotational volume accessible to the NC in its bound state (figure 5*c*), (4) *d*), (5) *e*) and (6) the dissociation constant *f*); *K*_{a} is defined in equation (2.6), which encodes the exact relation between these observables and the binding avidity.

From the main panels in figure 5*a*, it is evident that, unlike for flat substrates, *P*(*n*_{b}) for a cell membrane shows a peak at higher multivalency with a broader distribution. Both the flat substrates and the cell membranes show a distinct annulus-like pattern with the annulus size consistent with the respective multivalency distributions; this is shown in the insets to figure 5*a*. Since *P*(*t*_{1},*t*_{2}) is related to the lateral organization of the bound receptors, it is a determinant of the entropic contribution due to the receptors. This contribution is depicted in figure 5*b*: here *n*_{b}. We also find that the values of

In a similar manner, the rotational entropy may be estimated from the rotational volume Δ*ϕ*Δ*θ*Δ*ψ*, where Δ denotes the standard deviation in the three Euler angles *ϕ*, *θ* and *ψ* characterizing NC orientation about its centre of mass (see §2). The accessible rotational volume as a function of the receptor density is given in figure 5*c*, and we find that its estimate for both systems (cell membrane and flat substrate) are not sensitive to variations in *N*_{ant}.

We have also computed *d*, the value of *N*_{ant} which, as explained before, is due to the high affinity of the NC for the cell membrane even at low receptor densities. On the other hand, *N*_{ant}, which is reflective of a monotonic increase in the multivalency of the NC that is solely governed by the receptor expression levels (figure 5*a*).

Figure 5*e* shows the depth of the PMF *N*_{ant}, while its response in the case of cell membranes is non-monotonic. At low and intermediate surface densities, as is shown for membranes expressing 200–2000 ICAM1 molecules μ*m*^{−2}, *N*_{ant}, with the depth being maximum for *N*_{ant}=2000, and increases to a higher value when *N*_{ant}=4000. This behaviour indicates that a functionalized NC does not show enhancement in binding free energy if the target receptor expression exceeds a critical threshold value.

The above observations validate our hypothesis that the avidity of a functionalized NC is a result of the complex interplay between the various energetic and entropic contributions in the system. This is further exemplified in figure 5*f* where we depict the effect of *N*_{ant} and *κ* on the dissociation constant *K*_{d}, computed from the association constant *K*_{a} as *K*_{d}=(*K*_{a})^{−1}; for a definition of *K*_{a}, see equation (2.6). The values are reported as ratios normalized by *K*^{†}_{d}, which is the value of *K*_{d} for a system with *κ*=40 *k*_{B}*T*, *m*^{−2}, representative values for regular EC in lung. It is evident from figure 5*f* that the computed values of *K*_{d} are strongly influenced by *κ* and ICAM1 expression.

### 3.6 Predicting the endothelial targeting of ICAM1-coated nanocarrier in mouse

Figure 6*a*(i) shows the *in vivo* endothelial targeting of ICAM1-targeted NCs, with *N*_{ab}=41, 100 and 162/NC, in the lung, heart, kidney, liver and spleen of mouse [12]. These results show that a high level of targeting is achieved in the lung as evidenced by the increase in *η*_{exp} with increasing *N*_{ab}. On the contrary, the uptake in liver and spleen is non-targeted since *η*_{exp} is not sensitive to increase in *N*_{ab}. In panels (ii)–(v) of figure 6*a*, we present the corresponding computational predictions using four different model scenarios (see equation (2.7)) to quantify tissue targeting: (1) to flat substrates, (2) to endothelial cell membranes, (3) to endothelial cell membranes as well as accounting for capture by resting macrophages and (4) to endothelial cell membranes and accounting for capture by activated macrophages. To test the performance of our model, the corresponding Pearson's correlation coefficients [38] (denoted by *r*^{2}) from the comparison of the model results with experiments across all organs are given in the centre of each panel. Additionally, in order to evaluate the performance of our model in predicting the targeted contributions, we have computed the correlation coefficient for the lung alone, as this organ consistently shows the effect due to targeting, and the corresponding *r*^{2} values are shown in the top left corner of the panels in figure 6*a*. Our results show that the targeting behaviour is well captured by our model for NC binding to endothelial cell membranes (*r*^{2}∼0.88, panel (iii)) compared with the commonly used flat substrate models (*r*^{2}∼0.57, panel (ii)). This quantitatively verifies the hypothesis that the mechanical properties of the endothelial cell membrane is an important contributor to NC targeting. The inclusion of additional contributions due to the macrophages does not alter the targeting behaviour in lung (panels (iv) and (v)) as the targeting component of the avidity is the dominant term for mouse lung tissues.

The inclusion of capture by other cells (e.g. macrophages or monocytes) actually worsens the prediction accuracy of NC targeting levels in all organs. The correlation coefficient computed across all organs is highest (*r*^{2}=0.78) for the EC membrane model (i.e. without contribution from other cells). This implies that the non-specific uptake enhanced by EC capture is the dominant factor in differentiating the tissue targeting in the organs considered here. Hence, based on the statistical metrics we have presented in figure 6, we conclude that the model including the contributions from the endothelial cell membrane alone is the optimal choice in predicting NC targeting in mice tissues, and we use this model to make comparisons with other (additional) experiments (see below). We note that based on the above results, a very useful approximation for *K*_{EC}. The statistical performance of the various model scenarios including statistical significance of the model predictive power is summarized in table 3. The statistical procedures are explained in §2.6.

After having tested our model predictions against *in vivo* results for tissue targeting in mice, we now employ the model to make predictions for the tissue targeting of ICAM1-coated NCs in scenarios reported by several other studies [40,53–56]. The experimental values of *η*_{exp} are shown in figure 6*b* and, as noted earlier, the data show a large spread for spleen and liver. With the model parameters relevant for the tissue and NC in these experiments (see §2.7), we predict *η*_{sim}, using equations (2.6) and (2.7), and the results are compared with *η*_{exp} in figure 6*c*. Our model predictions with zero-fitted parameters show good agreement with experimental results (*r*^{2}=0.8), and points to the fact that our computational framework can give reliable and robust estimates for the tissue targeting of functionalized NC. The degree of agreement between our model predictions and six different sets of experimental results (figure 6) marks an important advance in the rational design of functionalized NC, with our proposed model framework providing a biophysical route for the optimization of functionalized NCs.

### 3.7 Predictions of tissue targeting of ICAM1-coated nanocarriers in human organs

After having tested and shown the predictive capabilities of the computational framework using data from several tissue targeting experiments in mouse, we now proceed to make predictions for NC targeting in human organs. We treat the endothelial cells in the liver, lung, heart, kidney and spleen of humans to be mechanotypically (i.e. the values of the mechanical variables of the membrane) similar to the corresponding organs in mouse. The protein expression levels in the various organs are determined, as before, using data from high-throughput mRNA and mass spectrometry measurements in human tissue reported in the BIOGPS and PAXDb databases. We note that this estimate is based on the caveat that cellular regulation of the target protein expression and accessibility are often physiology- and pathology-dependent and are often not known, even though the corresponding data in mice showed excellent correlation between gene expression and target expression (figure 3). The density of ICAM1 receptors in the five organs of interest are shown in figure 7*a* and we find the expression levels to be the largest in the spleen (approx. 2500 ICAM1 μ*m*^{−2}), while lung shows lower expression compared with that in mouse. The predictions for the targeting of anti-ICAM1-coated NCs, with *N*_{ab}=41,100 162/NC, are shown in figure 7*b*. Our results for ICAM1-targeted carriers indicate that the NC tissue targeting in human lung shows characteristic targeting behaviour but its sensitivity is less pronounced compared with that seen in mouse (figure 6*a*). By contrast, NC tissue uptake in the spleen is much larger compared with lung, where a distinctive signature of targeting behaviour is also predicted by our model. In other organs, the evidence for any targeting is absent in our model predictions. Owing to the large discrepancy between PAXDb and BIOGPS data for heart and liver (figure 7*a*), we have also analysed the sensitivity of the tissue targeting levels by varying the expression levels of *N*_{ant} compared with that obtained using expression levels in mouse. We find that the lung and spleen shows a pronounced sensitivity to variations in the receptor expression levels compared with other organs.

## 4. Conclusion

Rational approaches in the design of functionalized NCs for targeted delivery can be immensely beneficial in optimizing the affinity of the NC in a tissue and species specific manner. We have presented a molecularly guided and biophysically based model framework to predict live-cell/tissue targeting of functionalized NCs across multiple organs and species. This framework is based upon a *zero-fit biophysically based multiscale model* to compute the binding avidity for an NC binding to a cell surface. In our model, we distinguish the various target endothelial cells in terms of their *mechanotype* (i.e. the mechanical properties of the cell membrane) and *phenotype* (i.e. the expression levels of the target protein). We also account for the contributions from *non-specific mechanisms in the tissue*. The model presented here combines the model-predicted estimates of the partitioning coefficients for affinity–ligand functionalized NCs in the target tissue and the experimentally determined partition coefficient for a non-targeted carrier of the same size and type, hence taking both the targeted and the non-targeted contributions to tissue targeting.

Our results emphasize the fact that the mechanotype and phenotype of the target cell are key parameters that can significantly influence NC binding, and hence should be integral to the design of functionalized NCs. We have used our computational framework to predict the tissue targeting levels of 100 nm anti-ICAM1 functionalized NCs in the lung, heart, liver, kidney and spleen of mouse and compared our findings to *in vivo* experiments, where available. Our results show that the targeting behaviour of anti-ICAM1 functionalized NCs in the mouse lung can be well captured only if the mechanical properties of the endothelial cell membrane and entropic effects coupled to multivalent binding are both explicitly taken into account. Predictions for other organs, which show characteristic non-targeted behaviour in *in vivo* experiments, also depend critically on target cell mechanotype and phenotype. Furthermore, we have also tested the performance of our model in predicting tissue targeting levels in the various scenarios reported in five different experiments available in the literature, and as demonstrated by the statistical metrics of the comparison, our predictions agree very well with the experimental findings. We have also demonstrated how our model predictions can be tailored to other organisms and organs by computing the targeting of ICAM1-coated NCs in conditions that mimic human tissues. It should be noted that in terms of parameter determination, there is a larger variation between the physically based antibody-tracing method and other indirect methods of ICAM1-expression determination in spleen and liver than in the lungs. This suggests that the relative contribution of ICAM1 in the vascular lumen relative to the other compartments is organ specific. Hence, experiments directly mapping ICAM1 accessibility rather than indirectly determining gene or protein expression may be even more important for better characterization of the tissue microenvironment.

One of the limitations of our current model is that it only explores steady-state behaviour at timescales when internalization mediated by ICAM1 has not set in. We also note that our predictions rely on the accurate quantification of the ‘non-targeted’ capture (i.e. not directly mediated by ICAM1 binding) contributions through experiments since predictively accounting for such contributions is beyond the scope of the current model. Non-specific factors such as clearance mechanisms can be different for targeted and non-targeted uptake. The different cell types (including and beyond the endothelial cells and macrophages) in the tissue microenvironment can sustain context-dependent uptake mechanisms such as ICAM1- and FC receptor-mediated endocytosis, or macropinocytosis. Other physiological and haemodynamic factors such as vascular hydrodynamics and NC margination [16,64], glycocalyx [13,20], cytoskeletal dynamics and NC internalization [36,37] are also physiology- and pathology-dependent. These considerations warrant more work and will be addressed in the future through which we expect a significant improvement in *r*^{2} between model and experiment. Future extensions of the next-generation PK model presented here will also focus on predictively modelling the targeting of other epitopes in the vasculature (such as PECAMs), modelling targeting of epitopes in epithelial tissues and predictive modelling of non-targeted contributions discussed above.

## Ethics

*In vivo* studies in mouse were carried out in accordance with the Guide for the Care and Use of Laboratory Animals as adopted and promulgated by the US National Institutes of Health.

## Data accessibility

Sections S1–S4 contain detailed description of the computational model and derivation of the association constant and standardized uptake values. Movies M1–M7 show the binding of an NC to the seven different cell membranes shown in figure 4. These are available online under Supplementary Information. Access to raw data and codes can be requested through email to the corresponding author.

## Authors' contributions

N.R. and R.R. designed research, performed simulations and analysed data. R.W.T. implemented and obtained the results for the mean field cytoskeletal model. V.M. provided the in vivo experimental data. N.R., P.S.A., D.M.E., V.M. and R.R. discussed the results and wrote the paper.

## Competing interests

The authors declare no competing interests.

## Funding

Efforts to develop the membrane model for NC cell adhesion and the model for cytoskeleton were supported in part by grant nos. NSF: DMR-1120901, CBET-1236514. The free energy analysis for NC adhesion was supported in part by NIH: U01EB016027, and the pharmacodynamic model development was supported in part by NIH: 1R01EB006818-05. The *in vivo* experiments were supported in part by grant nos. NIH: HL125462 and HL087936. Computational resources were provided in part by the grant no. MCB060006 from XSEDE.

## Acknowledgements

We acknowledge helpful discussions from the members of the Radhakrishnan, Ayyaswamy and Eckmann laboratories.

- Received April 13, 2016.
- Accepted May 24, 2016.

- © 2016 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.