Originally published as Biophys J. BioFAST on November 3, 2006.
doi:10.1529/biophysj.106.095521
Biophysical Journal 92:831-846 (2007)
© 2007 The Biophysical Society
Understanding Intracellular Transport Processes Pertinent to Synthetic Gene Delivery via Stochastic Simulations and Sensitivity Analyses
Anh-Tuan Dinh,
Chinmay Pangarkar,
Theo Theofanous and
Samir Mitragotri
Department of Chemical Engineering, University of California, Santa Barbara, California
Correspondence: Address reprint requests to Prof. Samir Mitragotri, Dept. of Chemical Engineering, University of California, Santa Barbara, CA 93106. Tel.: 805-893-7532; E-mail: samir{at}engineering.ucsb.edu.
 |
ABSTRACT
|
|---|
A major challenge in synthetic gene delivery is to quantitatively predict the optimal design of polymer-based gene carriers (polyplexes). Here, we report a consistent, integrated, and fundamentally grounded computational methodology to address this challenge. This is achieved by accurately representing the spatio-temporal dynamics of intracellular structures and by describing the interactions between gene carriers and cellular components at a discrete, nanoscale level. This enables the applications of systems tools such as optimization and sensitivity analysis to search for the best combination of systems parameters. We validate the approach using DNA delivery by polyethylenimine as an example. We show that the cell topology (e.g., size, circularity, and dimensionality) strongly influences the spatiotemporal distribution of gene carriers, and consequently, their optimal intracellular pathways. The model shows that there exists an upper limit on polyplexes' intracellular delivery efficiency due to their inability to protect DNA until nuclear entry. The model predicts that even for optimally designed polyethylenimine vectors, only 
1% of total DNA is delivered to the nucleus. Based on comparison with gene delivery by viruses, the model suggests possible strategies to significantly improve transfection efficiencies of synthetic gene vectors.
|
INTRODUCTION
|
|---|
Polymer-based systems are being extensively studied as carriers for gene therapy (1
), and have also been used in clinical trials (2). Despite the obvious advantages of safety and malleability over viral vectors, polymer-based vectors (polyplexes) have not been very successful at the clinical level owing to their poor delivery efficiency. Even the best polyplexes are 1000-fold less effective than typical viral vectors such as adenoviruses. This poor efficiency stems from the numerous intracellular barriers that retain or destroy a majority of the gene dose before it can reach the host nucleus (3). A central challenge in the field is identification and synthesis of a polymer that can enhance translocation of the polyplex across these barriers. To this end, a large number (literally thousands) of polymers have been evaluated—poly-imines (4), dendrimers (4), polyamino esters (5), chitosans (6), and cyclodextrins (7), to name a few. However, only a few of these have been found to be significantly better than polyethylenimine (PEI25kDa), the accepted standard in polymer-based gene delivery. Further, it is not clear whether these polymers, many of which have been developed and optimized for use with cultured cells, will be effective in clinical applications.
The low transfection efficiency of synthetic vectors leads to a key question: Is there an inherent limitation to polyplexes? Or have we not found the magic polymer yet? Two challenges need to be addressed before this question can be answered. The first challenge, biological in nature, involves developing the design criteria of polyplexes that will lead to maximum efficiency. The second challenge, chemical in nature, involves design and synthesis of polymers that will form such polyplexes.
In this article, we exclusively address the first challenge. There are two major issues that make this task difficult. The first issue stems from fragmented understanding of intracellular trafficking of polyplexes, mainly due to the strong dichotomy in the experimental approaches used. One approach is based on macroscopic measurements of the delivery efficiency using in vitro transfection assays. The other approach is based on measurements at the single-cell single-particle level. Macroscopic experiments do not provide in-depth mechanistic understanding of the events leading to the final effect. Single particle data, on the other hand, aim to quantify individual transport steps (8,9), many of which, such as escape from endosomes, remain tremendously difficult to visualize due to the rare and random nature of these events (8). The link between these two distinct scales (i.e., macroscopic and microscopic) is still missing. The second outstanding issue is the limited knowledge of the structure-activity relationships that link vector's physico-chemistry to the efficiency of individual trafficking steps. This limitation confounds systematic experimental optimization of vector properties.
We address the problem of understanding and optimizing intracellular transport by developing a detailed mathematical model of the gene delivery pathway. Previous models of gene delivery have essentially followed a pharmacokinetic approach (10–13), which treats the cell as a well-mixed compartment. The major shortcoming of kinetic models is that they approximate all spatial and transport processes by kinetic equations. The kinetic rate constants have to be obtained by fitting the model to experimental data. This makes it difficult to extrapolate the results of the model predictions. Most importantly, kinetic models, due to their simple approximations of intracellular transport, fail to take advantage of the wealth of data made available by single-particle tracking experiments.
We realize that much of the complexity of the vector-cell system arises from the spatio-temporal variation in the rates of various intracellular processes and can be captured only by considering a spatial view of the cell. To address this challenge, we develop an advanced stochastic simulation framework that effectively describes the dynamic and spatial nature of intracellular trafficking of nanoscale carriers. By providing a realistic representation of topology and dynamic organization of the cell interior, the modeling framework serves as a three-dimensional, "live," computer-constructed cellular map for navigation of nanoscale carriers. The model contains information from the shortest to the longest relevant length and timescales, and wherever possible, is rigorously validated against experimental data. Wherever mechanistic or quantitative information in the literature is incomplete, we perform systematic experiments to measure the required processes. We use the model to calculate the delivery efficiency—the probability that a polyplex, upon internalization, successfully delivers DNA to the host nucleus, and to predict the optimal intracellular itinerary. We show that the optimal pathway is controlled by several cell-specific properties such as cell morphology, endocytic trafficking, and microtubule-dependent transport.
|
EXPERIMENTAL METHODS
|
|---|
Polyplexes
Polyethylenimine (25 kDa) was obtained from Sigma (St. Louis, MI). All fluorescent dyes (Oregon green, TMR-Dextran, 70 kDa) were obtained from Molecular Probes (Eugene, OR). Amine groups on polyethylenimine (25 kDa) were labeled using the succinimidyl ester of Oregon green as per the manufacturer's protocol. The complexes were always prepared at a Nitrogen/Phosphorus ratio of 9, in 150 mM NaCl solution. Complexes were used between 30 and 45 min post-synthesis. They were verified to transfect skin fibroblasts.
Cell culture and transfection
Human Skin Normal fibroblasts (ATCC, TE.353.Sk, Manassas, VA) were cultured as per the standard ATCC protocols. Briefly, cells were maintained in Dulbecco's modified eagle's medium (ATCC), which had been supplemented with 4 mM L-glutamine and 1.5 g/L sodium bicarbonate. Cells were grown at 37°C and 5% CO2; and maintained in the log phase by consistent subculture. On the day before the experiment, cells were seeded onto a glass bottom petri dish at a density of 7000–10,000 cells/cm2. The next day, complexes were added to serum-free medium at a concentration of 2 µg DNA/ml. The cells were incubated with this medium for 15–30 min at 37°C. Post-incubation, the unbound particles were washed away 5x by phosphate-buffered saline and the cells were replenished with fresh serum containing medium, which was supplemented with 25 mM HEPES to maintain carbonate-bicarbonate balance in the absence of CO2.
Fluorescence microscopy
The cells were placed on the stage of a Zeiss (Thornwood, NY) Axiovert-25 inverted microscope, which was fitted with a Bioptechs Delta T temperature controller (Bioptechs, Butler, PA). The cells were observed using an oil immersion 100x objective. Oregon green was excited and observed using HQ480/40x and HQ 555/40M FITC compatible filters from Chroma (Rockingham, VT). For every cell, bright-field and fluorescence images were acquired using a cooled CCD-camera (CoolSNAPHQ, Roper Scientific, Duluth, GA) which was controlled by Metamorph imaging software (Universal Imaging, Downington, PA). The images were acquired for 1–2 min at 1 frame per second.
Immunofluorescence
Lysosomes in the cells were immunolabeled using standard techniques. Briefly, the cells were fixed for 30 min using freshly prepared 2% paraformaldehyde. This was followed by permeabilization with 0.01% Tween-20 in phosphate-buffered saline and blocking with 2% BSA for 2 h. The cells were then labeled with mouse monoclonal anti-Lamp1 (primary antibody, BD Biosciences) and subsequently with Rhodamine red-X-Goat anti-mouse IgG (secondary antibody, R6393, Molecular Probes). The cells were then mounted using VectaShield (Vector Labs, Burlingame, CA). Rhodamine-X labeled antibodies were observed using HQ 545/30X and HQ 640/50M filters from Chroma.
Image processing and analysis
All image acquisition and processing steps were performed within the Metamorph software. Phase contrast images of cells were used to extract coordinates of the cell membrane and nuclear membrane. The raw fluorescence images show bright punctuate structures corresponding to diffraction-limited images of fluorescent vectors. Images were processed with a custom FFT-based protocol to filter out very long-range and very short-range correlations. This was followed by a series of Gaussian blurring, convolution, and unsharp masking steps until the resultant image could be thresholded such that most of the bright structures could be individually measured. All particle coordinates were exported to MS Excel (Microsoft, Redmond, WA). All further calculations involving cell shape, vector distribution, etc., were performed in MatLab (The MathWorks, Natick, MA). For single particle tracking, stacks of fluorescent images of cells were filtered to remove high frequency noise. Other filtering steps were avoided to preserve the intensity structure of the raw image. Trajectories of all the particles were individually verified by replaying the video with the trajectory superimposed on the image. The coordinates of a vector in all the frames were exported to MS Excel, and were used for all further calculations. Colocalization of Oregon green labeled PEI25-DNA vectors and anti-Lamp-1 labeled lysosomes was measured by overlaying the Red and Green channel images. The images were self-normalized using standard histogram-stretching methods. Colocalization was considered significant if the Red and Green structures had >85% of their combined pixels in common.
|
STOCHASTIC SIMULATIONS
|
|---|
The objective of stochastic simulations is to bring-out the spatial and physical aspects of intracellular trafficking of polyplexes by providing a realistic representation of cell geometry and intracellular organization and a discrete and mechanistic description of transport processes. As in single-particle tracking experiments, the proposed model follows the trajectories of polyplexes inside cells. The computational domain is defined by the cell geometry, which is reconstructed from experimental visualization of cells (see below). We distinguish between in vitro and in vivo experiments. Many cells, especially those considered in this analysis (human skin fibroblasts) spread on the substrate under in vitro conditions, and the cell thickness is often much smaller than other dimensions due to adhesion to the surface. Therefore, it is sufficient to formulate the simulations in two dimensions for in vitro applications (Fig. 1 a). In vivo applications, however, require considerations of three-dimensional geometry (Fig. 1, b–d).
View larger version (91K):
[in this window]
[in a new window]
|
FIGURE 1 Computer representation of human skin fibroblasts (a) in vitro, where the cell is represented as a two-dimensional object and (b) in vivo, where the cell is represented by a three-dimensional structure. Green and red objects represent endosomes and lysosomes respectively. Thin yellow lines represent microtubules (MTs). The microtubule-organizing center (MTOC) is randomly located in the vicinity of the nucleus. For clarity, we use fewer endosomes, lysosomes, and MTs to render the reconstructed cells than in simulations. Two regions of the three-dimensional cell are magnified and shown in panels c and d.
|
|
The positions and the biophysical states of polyplexes are updated using a stochastic algorithm. Fundamentally, all stochastic events simulated here can be classified into two categories: 1), transport events; and 2), reaction events. Transport events (diffusion or motor-driven movements) are captured by equations of motion. Reaction events, involving association, dissociation or rupturing of entities, are approximated as first-order Markov processes. Definitions of variables and symbols are summarized in Table 1.
Mathematical description of cellular environments
Cell geometry
For in vitro applications, we reconstruct cell geometries from phase contrast images of cultured human fibroblasts. Viewing from the top, the cell interior is divided into two regions, cytoplasmic and supranuclear, separated by the nuclear boundary (see Fig. 1 a). We use an image analysis program (Metamorph, Universal Imaging) to extract the coordinates of the nuclear boundary and the cell boundary. The high-resolution images captured at 100x magnification allow for precise determination of the shape of cell and nucleus. For in vivo applications, we reconstruct the cells from histological and electron microscopy images of fibroblasts embedded in tissues. The procedure for in vivo cells is more complicated and described in detail in the Supplementary Material. The cell geometry is represented by two surfaces describing the cell membrane and the nuclear envelope (see Fig. 1 b). Our reconstructed three-dimensional cells are in good agreement with past descriptions of fibroblasts under in vivo conditions. All cell geometries remain unchanged during simulations.
Microtubules
Microtubules (MTs) exhibit dynamic instability, i.e., stochastic switching between prolonged phases of assembly and rapid phases of disassembly. The frequency of transition from slow assembly to rapid disassembly (referred to as "catastrophe") determines the stability of the MT population (24). In this model, MTs are modeled as straight lines, which grow radially in random directions from the microtubule-organizing center (MTOC) toward the cell periphery. MTOC is placed randomly in the vicinity of the nucleus. To describe microtubule polymerization dynamics, we use stochastic equations that follow the individual MT length history. A microtubule switches intermittently between three phases (growth, shrinkage, and collapse) as follows (14):
 | (1) |
where q and m are random numbers between 0 and 1. The values kgrow, kshrink, and kcat are the reaction rate constants for MT growth, shrinkage, and collapse, respectively.
Mathematical description of intracellular trafficking of polyplexes
We introduce two state variables to characterize intracellular trafficking of polyplexes, namely, biophysical state S, which indicates the biological compartmentalization and transport state s, which represents the type of movements of polyplexes.
Biophysical states
After binding to the cell membrane, the intracellular pathway of a polyplexes can be represented by five distinct biophysical states S: bound to cell membrane (M); inside endosomes (E); inside lysosomes (L); inside cytoplasm (C); and unpacked plasmids (P) (Fig. 2). The states correspond to the stages of the gene delivery pathway (see Fig. 2 legend and Supplementary Material, section 1.3).
View larger version (18K):
[in this window]
[in a new window]
|
FIGURE 2 Gene delivery pathway of polyplexes. Polyplexes are internalized by endocytosis. They are immediately delivered to endosomes and subsequently to lysosomes. To gain entry to the nucleus, the polyplexes must escape from endosomes or lysosomes. Once inside the cytoplasm, the polyplexes unpack to release the DNA for successful nuclear entry. Ultimately, the exogenous DNA is transcribed, and gene products are synthesized. The biophysical states S and the rates of transition between them are shown in the figure.
|
|
Transitions from one state to another depend on various biological and physical factors, such as location, interactions with cellular organelles, and physiochemical properties of polyplexes. In the present model, we approximate the transition between two biophysical states of the vectors as a first-order reaction. For instance, we use kescape (i.e., the rate of escape from endocytic vesicles) to characterize the transition from S = E to S = C (Fig. 2) The probability that a transition from state E to state C occurs in the interval [t, t+
t] is given by
 | (2) |
We then compare a random number uniformly distributed between 0 and 1 to this probability to determine whether or not the event occurs in the time interval t (see Supplementary Material, section 1, for more information). Other transitions are predicted in a similar manner.
Transport states
Each biophysical state has a distinct transport pattern. Membrane-bound and cytoplasmic polyplexes move only via diffusion. On the other hand, polyplexes inside endosomes switch intermittently between diffusion and directional transport. To account for these distinct transport modes, we introduce the concept of transport state s. Each biophysical state is characterized by a set of transport states and a diagram depicting the transitions between these states. At each point in time, a polyplex occupies one of the movement states. For instance, transport of polyplexes inside endosomes is represented by three transport states, s = ±1, –1, and 0 (+1, plus-end directional transport; –1, minus-end directional transport; and 0, diffusion). We define 
and 
as velocities of directional movements toward the plus-end and the minus-end of microtubules. The values 
and 
are first-order rate constants characterizing particle's transition from s = 0 to s = +1 and s = –1, respectively. The values 
and 
are the rate constants of corresponding reverse processes. Thus, the cytoplasmic transport pattern of endosomes is represented by seven parameters, 
(free diffusion coefficient), 
, 
, 
, 
, and 
. These transport coefficients are independent of location inside the cytoplasmic region. Quantitative measurements of movements of PEI-DNA complexes in human fibroblast support this assumption. Table 2 summarizes the transition diagram and the associated transport parameters for all five biophysical states.
The position of a polyplex occupying a diffusive transport state is updated using two- or three-dimensional random walk. For MT-bound entities, the equations of motion are applied,
 | (3) |
where 
is the position of polyplex i and 
characterizes the microtubule j that the polyplex i is associated with.
To simulate transitions between biophysical states, we employ the same principle introduced above (see Eq. 2). However, the occurrence of several transition events requires specific physical conditions. For instance, binding of a diffusing endosome to MT only occurs if the endosome is within a radius re from a neighboring MT, and the probability that it will bind to and travel toward the plus- and minus-ends of the MT in an interval t are 
and 
, respectively. Similarly, termination of microtubule-dependent runs depend not only on the rates of unbinding, and , but also on the physical layout of the microtubule tracks. The endosome will fall off the MT, that is, become diffusive in the next time step, once reaching the ends of the MT or hitting the nuclear boundary or the cell boundary, which act like solid walls.
The most interesting case is the transition between the clustered state and the free state of polyplexes inside lysosomes. Lysosomes tend to form clusters in the perinuclear region (Fig. 4 a, inset). When a lysosome-carrying polyplex aggregates with other lysosomes, it exhibits restricted diffusion. When it is free, it can travel bidirectionally along microtubules just like endosomes. The declustering rate, 
, is presumably space-dependent, while the rate of transition from the free to the clustered state is a function of the normalized clustering rate and the local concentration of lysosomes, CL(r), which is obtained a priori by analysis of lysosome images after immunolabeling.
View larger version (16K):
[in this window]
[in a new window]
|
FIGURE 4 Intracellular spatial distribution of lysosomes and transfer of PEI-DNA from endosomes to lysosomes. (a) Fluorescent micrograph of lysosomes labeled by incubation of human dermal fibroblasts with fluorescent dextran for 2 h followed by a 22-h chase in culture medium. The nuclear and cell membrane are traced for clarity. Inset shows typical clusters of perinuclear lysosomes. (b) Spatial distribution of lysosomes represented by p(r), the probability of finding a lysosome at a distance r from the nucleus. (c) Gradual delivery of PEI25-DNA vectors from endosomes to lysosomes. The ordinate is the fraction of total intracellular PEI25-DNA vectors that colocalize with LAMP1 positive structures (lysosomes) at different times post-transfection. Lysosomal delivery is seen to be substantially complete at 24 h post-transfection. The shaded line is a first-order kinetic fit on the data yielding ktransfer = 1/500 min–1.
|
|
Simulation scheme and data analysis
Here, we describe the simulation scheme for in vitro applications. The simulation scheme for in vivo applications is similar. For each simulation, an arbitrary cell configuration is chosen from the library of different cell geometries, which defines the computational domain. During one time-step, the cellular environments and the states/positions of polyplexes are updated according to the following scheme:- Update cellular environments: cell geometry (if desired), and the state (growth, shrinkage, or collapse) and the length of microtubules.
- Updates of polyplexes: biophysical states, transport states, and polyplexes' locations.
To guarantee accuracy and stability of the numerical scheme, we use a small time step, t = 0.2 s. The time step is much smaller than the governing timescales of the transport processes. Varying t between 0.1 s and 2 s does not significantly influence the model predictions. For each parameter set, we sample 10,000 trajectories from 200 cells. The delivery efficiency is calculated based on the number of DNA molecules that gain entry to the nucleus at a final time Tf, when we terminate the simulations. The simulation software is available by request.
|
PARAMETER ESTIMATION
|
|---|
The model contains a large number of parameters, which describe different aspects of the problem such as cell geometry, microtubules, trafficking of endocytic vesicles, trafficking of vectors in the endocytic pathway, and trafficking of vectors after endosomal/lysosomal escape. We measured parameters describing the trafficking of vectors present within endocytic vesicles using skin-fibroblasts as the model cell-line and PEI25kDa-DNA as the vector. However, for the steps after endosomal escape, experimental measurements are not possible at the single particle level. In such cases, data from bulk measurements was used. For example, we adapted the rate of plasmid degradation in the cytoplasm from Lechardeur et al. (15), who measured degradation of injected plasmids using fluorescence in situ hybridization.
Microtubule organization and dynamics
Based on immunolabeled images of fibroblasts, we estimate that the number of microtubules in the plane of movements under in vitro conditions to be 1000–1500. To our knowledge, microtubule polymerization dynamics in human skin fibroblasts have not been quantified yet. Accordingly, we postulate that the parameters that characterize microtubule dynamic instability in human fibroblasts are similar to those of CHO fibroblasts. Komarova et al. (16) investigated the lifecycle of MTs in CHO cells and reported a growth rate (vgrow) of 0.25–0.35 µm/s, a shortening rate (vshrink) of 0.3–0.4 µm/s, and a catastrophe frequency (fcat) of 0.005 s–1. Also, the mean growth time tgrow is 40–60 s and the mean shrinkage time tshrink is 5–15 s. Variations of vgrow, vshrink, tgrow, and tshrink within the estimated intervals do not significantly affect the predicted particle distributions. The most sensitive parameter is fcat, which determines the stability of the microtubule tracks.
Transport properties of endosomes
The cytoplasmic transport properties of endosomes are represented by seven parameters, D0,E (free diffusion coefficient), 
, 
, , , , and . These parameters are not available in the literature and are therefore measured from our own experiments (see Supplementary Material, Section 1.2.1).
- Diffusivity, D0,E. Trajectories of 50 fluorescently labeled, almost immobile vectors were selected from experimentally measured movements of PEI25kDa-DNA vectors (assuming that vectors are present within endosomes). These trajectories were used to calculate the mean-squared displacement (MSD) as a function of the timescale t. Based on the MSD, we estimate D0,E to be 0.0005–0.001 µm2/s. Similar values have been reported previously in the literature. Immunolabeling of endosomes (using anti-rab5, red) showed that all observable green fluorescence of PEI-DNA vectors colocalizes with anti-rab5 at the time of measurements (1 h post-transfection). This confirms that, within experimental errors, the imaged vectors were indeed present within endosomes.
- Rate of unbinding,
. The rate of unbinding is the inverse of the mean duration of "spurts" or "runs" in a certain direction. Fig. 3 a shows the probability distribution of duration of a run in the plus and minus directions, calculated from 250 active trajectories. The exponential nature of this distribution also supports the assumption of a first-order unbinding process. The run lengths in either direction show identical distributions and thus yield 
.
- Rate of binding,
. At equilibrium, a certain fraction f of particles must exist in bound states 1 and –1, and 
. Here 
and 
; and for 
, we get 
. Thus, knowing 
, and measuring f from time-lapse videos of cells at different times, the rate of binding can be calculated. For PEI-DNA containing endosomes, the rate of binding calculated from videos over 20–120 min post-transfection is 1/8 s–1.
- Velocity distribution,
. Endosomes do not exhibit a constant velocity, but a distribution of velocities that ranges from <0.05 µm/s to 3 µm/s. Fig. 3 b shows the distribution of frame-to-frame speeds in both the plus and minus directions. It can be seen that the distributions are identical exponential distributions, with a mean velocity of 0.33 µm/s.
View larger version (8K):
[in this window]
[in a new window]
|
FIGURE 3 Transport properties of endosomes. (a) Distribution of the duration of runs of endosomes in the plus (solid circles) and minus (open circles) directions and (b) distribution of the frame-to-frame velocities in plus (solid circles) and minus (open circles) directions. The vertical coordinate indicates the frequency of occurrence of an event, i.e., the probability of observing a certain run-length (a) or a certain velocity (b).
|
|
Transport properties of lysosomes
Lysosomes form clusters and accumulate in the perinuclear region, as illustrated in Fig. 4, a and b. Lysosomes switch stochastically between free and clustered states, and the properties in each of these states are measured as follows:- Transport of free lysosomes. We use the same method as in case of endosomes to obtain
, 
, 
, and Df,L for free lysosomes. Lysosomes are labeled by fluorescent dextran (10 kDa), which has been chased into the lysosomes by incubating the cells for 24 h. We confirmed this experimentally by labeling endosomes with lysosome-associated-membrane-protein (lamp)-1. The analysis of the experimental data on lysosome transport yields = 0.33 s–1, = 1/20 s–1, = 0.56 µm/s, and Df,L = 10–4 µm2/s.
- Transport of clustered lysosomes. Clustered perinuclear lysosomes rarely exhibit directed motion. In the event that directed motion is seen, it is always followed by dissociation of the lysosomes from the cluster into the free phase. Thus, transport of a lysosome in the clustered phase is described by diffusion (Dc,L = 10–4 µm2/s).
- Transition between free and clustered states. To find the rate at which a free lysosome becomes part of a cluster and vice versa, we examined videos of lysosomal clusters. In every video, a cluster of 4–5 lysosomes was randomly chosen and tracked for a period of 2 min. The cluster is often highly dynamic, with lysosomes constantly leaving the cluster. At the same time, other lysosomes come to join the cluster. Very often, the departed lysosomes rejoin the original cluster. Averaging over 75 clusters, we obtain the rate of clustering, i.e., the rate of forming physical bonds between two lysosomes, to be kcluster = 0.5 s–1. Based on the global state of the cell (the number of free and the number of clustered lysosomes) and assuming equilibrium, we can find the equilibrium constant for the clustering-declustering processes,
= kcluster/kde-cluster. Here = 4, and hence the rate of declustering is kde-cluster = 0.125 s–1.
Transfer of polyplexes from endosomes to lysosomes
When polyplexes are inside endosomes and lysosomes, they display transport patterns according to the parameters reported above. To determine the rate of transfer from endosomes to lysosomes, ktr, we examine the colocalization of Oregon green labeled PEI-DNA (in endosomes) with rhodamine-immunolabeled lysosomes as a function of time, post-transfection. Quantitative measurements allow us to capture the increasing fraction of PEI-DNA present in lysosomes (Fig. 4 c). The data are averaged over each cell, and over multiple cells at each time point. The curve is modeled with a first-order reaction to yield ktr = 1/500 min–1. This is the global transfer rate. The local transfer rate will depend on the local concentration of lysosomes, cL.
Physiochemical properties of polyplexes after endosomal escape
The parameters which describe the processes after endosomal escape are adapted from the literature, including the rate of DNA degradation (15), the diffusivity of vector in the cytoplasm (17), the rate of unpacking (18), and the rate of binding of the vector to the nucleus (19) (see Supplementary Material, Table S4).
|
MODEL VALIDATION
|
|---|
The stochastic model spans several length and timescales, as discussed above. The model is therefore validated at three separate levels against: 1), exact solution of single vesicle transport on single filament at timescales of 100 s; 2), experimental measurements of mean-squared displacements of endosomes in a live cell at timescales of 100 s; and 3), whole-cell-level spatial distribution of PEI25kDa-DNA vectors over 24 h. Comparisons of model predictions to gene expression data are given in Supplementary Material, section 2.5.
Comparison with diffusion-reaction-advection models
To verify that the stochastic algorithm was implemented properly, we simulate movements of endosomes along a single, fixed filament. Transport of endosome obeys the transition map depicted earlier (see Table 2). The endosome is initially located at x = 0. The spatial movements of the endosomes in time are characterized by p(x,t), the probability that an endosome is located at position x at time t. This probability can be theoretically predicted by a system of one-dimensional diffusion-advection-reaction equations (20). Readers are referred to Supplementary Material in section 1.7 for more information on the diffusion-advection-reaction equations. In Fig. 5 a, we compare p(x,t) predicted by stochastic simulation (sampling of positions of 5000 particles) to the exact solution provided by diffusion-reaction-advection equations, for t = 50 s and 100 s. The good agreement between simulation and analytical solution confirms that the stochastic algorithm is capable of accurately predicting transport mediated by molecular motors.
View larger version (12K):
[in this window]
[in a new window]
|
FIGURE 5 Model validation. (a) Displacement probability p(x,t) of a particle engaging in motor-assisted transport along a filament at t = 25 s (solid line, exact solution; open circles, simulation) and t = 100 s (dashed line, exact solution; solid circles, simulation). (b) Average mean-squared displacement for stochastically generated trajectories (solid line) and the experimentally measured trajectories (triangle). A crossover from subballistic to diffusive regime occurs between t = 10–100 s.
|
|
Comparison with measured mean-square displacements of endosomes
We simulate movements of endosomes in two-dimensional cultured fibroblasts. The cell geometry is accurately reproduced from phase contrast images, and the MTOC is assigned a random location near the nucleus. The microtubules grow radially from the MTOC and exhibit dynamic instability, thus mimicking the real cells. A comparison between the predicted MSD and the measured MSD is shown in Fig. 5 b. At short times (t < 1 s), simulations show that vesicle motion is ballistic, as also seen in experiments (t
, 
2). However, at long times (t > 10 s), particle motion exhibits diffusive scaling (t1) and may be described by an effective diffusive coefficient, Deff . Kulkarni et al. reported similar observations (21). Over short timescales, vesicles experience limited runs and the mean velocity is nonzero and directional. However, at long times, vesicle motion is averaged over several random (in magnitude and direction) ballistic runs and pauses, and exhibits a diffusionlike behavior. The simulation captures the essential physics of microtubule-dependent transport.
We observe that Deff depends strongly on the density and the stability of microtubules. Decrease in the number of microtubules, NMT, and increase in the frequency of catastrophe, fcat, substantially reduce Deff and cause accumulation of endosomes in the periphery. The effects of microtubule organization on gene delivery have been discussed in detail elsewhere (22). It should be noted that the effects of microtubule density and stability on the overall delivery efficiency are not significant for the ranges of parameters under consideration.
Comparison with measured spatial distribution of PEI25kDa-DNA vectors—phases of perinuclear accumulation—cytoplasmic versus supranuclear vectors
The model was applied to predict the spatio-temporal distribution of PEI25kDa-DNA complexes that are present within endocytic vesicles. Prediction of the spatial distribution of vectors in dermal fibroblasts at 30 min, 4 h, and 11 h post-transfection (Fig. 6, lower panel) shows excellent qualitative agreement with experimental images (Fig. 6, upper panel). Looking from the top, the cell cytoplasm was divided into two regions—1), supranuclear; and 2), cytoplasmic—which are separated by the nuclear boundary. The supranuclear region is a thin layer of 2–4 µm in width directly above the nucleus. Polyplex distributions in these two regions were treated separately.
View larger version (59K):
[in this window]
[in a new window]
|
FIGURE 6 Comparison of model predictions (lower panel) of spatiotemporal distribution of PEI-DNA complexes to experimental data (upper panel). PEI-DNA complexes are labeled with Oregon green (green) and lysosomes are immunolabeled with rhodamine-X labeled secondary antibody (red). The yellow regions represent PEI-DNA vectors that colocalize with lysosomes. The same cells are reconstructed using simulations and presented in lower panel. The three cells correspond to (a, b) 30 min, (c, d) 4 h, and (e, f) 11 h post-transfection.
|
|
We also performed a quantitative comparison of model predictions with experimental data by calculating normalized polyplex concentration c(r,t) (the number of polyplexes per unit area at a distance r from the nuclear membrane at time t post-transfection). The predicted polyplex concentration is within ±10% of the experimental data (Fig. 7). This is encouraging since we do not use any fitting and all parameters are either measured or adapted from literature.
View larger version (14K):
[in this window]
[in a new window]
|
FIGURE 7 Validation of the model against experimental measurement of spatio-temporal distribution of PEI25kDa-vectors. (a) Normalized local concentration c(r,t) of vectors measured experimentally (symbols) and predicted by the model (lines) at 10 min (circles), 2 h (squares), and 24 h (triangles). The abscissa is r, the distance from the nucleus. The distributions are sampled from 20 cells. (b) The temporal evolution of vector concentration in the supranuclear region (squares) and the perinuclear region (circles). The perinuclear region is defined as a region of width 2 µm around the nucleus. The lines are the predicted distributions.
|
|
Initially, the polyplexes are dispersed uniformly in the cell due to random sites of endocytosis on the cell membrane (c(r,t = 0) = 1, Fig. 7 a) and accumulate near the nucleus with time. Perinuclear accumulation of polyplexes occurs in two distinct phases (Fig. 7 b, cperi(t)). During the first phase, the average concentration of complexes in the perinuclear region quickly doubles within 1–2 h. This rapid accumulation is due to facilitated diffusion of endosomes containing PEI-DNA complexes on MTs as discussed above. Consequently, the complexes gradually disperse along the MTs. Since MTs are organized in an asterlike fashion, a uniform distribution on MTs leads to a higher concentration of polyplexes in the perinuclear region (23). The timescale over which the first phase of accumulation occurs depends upon the size of the cell and the MT-based effective diffusivity. The second phase of perinuclear accumulation is much slower, with a timescale of 6–10 h, and arises from gradual transfer of complexes from endosomes to lysosomes (Fig. 4 c). Delivery to lysosomes immobilizes (and hence concentrates) the polyplexes in perinuclear clusters (Fig. 7 b, cperi(t)). Thus, perinuclear accumulation of complexes arises not from any active targeting, but from the transport processes inherent to endocytic trafficking. The model, based on fundamental interactions between vesicles and microtubules, is able to predict the spatio-temporal distribution of polyplexes, which evolve at much longer timescales.
In the supranuclear region, the density of MTs is low and most polyplexes do not exhibit MT-based transport. Supranuclear endosomal polyplexes diffuse slowly toward the perinuclear region, are transferred to lysosomes, and become immobilized in the perinuclear clusters. Effectively, polyplexes are transported out of the supranuclear region (Fig. 7 b, csupra(t)). The model accurately predicts the temporal evolution of vector concentration in the supranuclear region. The excellent agreement between the measured and calculated values of c(r,t), cperi(t), and csupra(t) is a strong indication of the model's capability to predict transport of polyplexes inside cells.
|
MODEL PREDICTIONS
|
|---|
The model was first used to examine the individual steps of the gene delivery pathway. The entire pathway was divided into two stages: 1), before endosomal/lysosomal escape, when the transport of vectors is determined by the endocytic vesicles; and 2), after endosomal escape, when the transport is determined by diffusion and stability of the vector. We first calculate the optimal location and time of endosomal escape that maximizes the probability of completing DNA delivery for in vitro cells. We then use this knowledge to predict the optimal pathways for in vitro and in vivo cells. Due to the complexity and multidimensionality of the parameter space, a global optimization of all parameters involved is not possible. Instead, we develop a more intuitive procedure to find the optimal set of parameters.
Optimal location of endosomal escape
After escape from endocytic vesicles, polyplexes must unpack to release the plasmid, which can then enter the nucleus. Diffusion coefficients of polyplexes and plasmids in the cytoplasm are small, 10–4 µm2/s and 10–3 µm2/s, respectively (17). While the polyplex is chemically stable, the plasmid is susceptible to cytoplasmic nucleases and has a short half-life (50–90 min) (15). Thus, only those plasmids that are close to the nucleus have a reasonable chance of success. To quantify this effect, we calculate the probability, 
(r), that a polyplex, having escaped from an endocytic vesicle at distance r from the nuclear boundary at time t, will successfully deliver DNA to the cell nucleus (Fig. 8 a) (see Supplementary Material for mathematical definition of (r)). The success probability (r) decreases as vectors escape away from the cell center, which is located at r = –5 µm. The value (r) is generally high for supranuclear escape (r < 0) and declines only slightly toward the nuclear boundary (r = 0 µm). On the contrary, for vectors that escape within the cytoplasmic region (r > 0), the probability decreases exponentially with r, and at r =1 µm, (r) is already reduced to mere 0.01 (that is, 1%). The prediction that success probability is low for vectors that escape far away from the nucleus is intuitive. However, predictions in Fig. 8 a provide significant new information. First, they provide a quantitative determination of how success probability decreases with distance. Second, they provide a direct means of comparing the contribution from supranuclear and cytoplasmic vectors. This will become important in discussing in vitro transfection experiments in which the supranuclear region often occupies a significant fraction of the cell (discussed later).
Optimal time of escape
A similar analysis was performed with respect to time of escape, te. Specifically, we calculated the probability 
(te) that a vector escaping at time te post-transfection reaches the nucleus within 24 h (see Supplementary Material for mathematical definition of (te)). (te) is essentially a spatially averaged value of (r) for all vectors released in the cytoplasm at time te. As seen in Fig. 8 b, (te) is very different for supranuclear and cytoplasmic vectors. For supranuclear vectors, (te) decreases monotonically with te as the vectors are gradually transported out of the supranuclear region due to endosome-lysosome fusion (Fig. 7 b, csupra). For cytoplasmic vectors, (te) exhibits an optimum in te at 12–13 h post-transfection. If endosomal escape occurs too early, when most of the cytoplasmic vectors are far away from the nucleus, the released vectors will degrade before reaching the nucleus. On the other hand, late escape also results in low delivery efficiency. First, the longer the vectors stay inside endocytic vesicles, the more likely they will unpack and degrade inside lysosomes. Second, the released vectors may not have sufficient time to find and to enter the nucleus before the simulation is terminated. The precise contribution of each effect to the existence and the actual value of optimal escape time depend on several parameters.
The data in Fig. 8 b provides several new insights into transport of gene vectors. First, it clearly establishes that the best strategy for supranuclear vectors is to escape as early as possible. On the other hand, there exists a clear optimum escape time for cytoplasmic vectors. But most importantly, the data reveals that due to the stark differences between the behavior of supranuclear and cytoplasmic vectors, the geometry of the cell must be taken into account while determining the optimal escape time. Different cell lines have different sizes of nucleus and nucleus area to cell area ratios. These factors must be taken into account while comparing data from different cell lines.
Global optimization of vector trafficking in vitro
We next combined the data from the previous two sections and determined an overall efficiency 
, the probability that a vector released anywhere within the cell at any time delivers DNA to the nucleus. We then sought to determine vector parameters that optimize . We first defined the possible ranges of the parameter values that polyplexes are expected to possess. The parameter space is highly multidimensional, and consequently full optimization is computationally demanding. Instead, we followed a more intuitive path to reach the global maximum. The delivery efficiency exhibits a simple monotonous relationship with many parameters, e.g., rates of internalization and DNA degradation. We maximized these parameters to reach a semioptimal state, around which we performed local optimization within a lower dimensional parameter space. We focused on a three-dimensional region of the parameter space represented by kescape (the rate of vector escape from endosomes or lysosomes), ktr (the rate of vector transfer from endosomes to lysosomes), and kunpack (the rate of vector unpacking, and consequently plasmid degradation either in lysosomes or in the cytoplasm). These parameters account for the delicate balance between endocytic trafficking, endosomal escape, and vector stability. All other vector parameters were held at their respective optima. It is important to note that the parameters under consideration are both vector-specific and cell-specific. For instance, our experiments show that the rate of transfer from endosomes to lysosomes in human fibroblasts, ktr, depends strongly on the nature of the internalized cargo: ktr for dextran, poly-lysine-DNA vectors, and PEI25kDa-DNA vectors is 1/45 min–1, 1/60 min–1, and 1/600 min–1, respectively. Molecular conjugation of gene vectors with proper ligands will significantly alter the value of ktr. The other parameters, e.g., kdegL, kescape, and kunpack, show the same dependency on the physicochemical properties of the vectors.
Under in vitro conditions (flat, adherent, two-dimensional geometry), optimization yielded relatively straightforward results, that is, vectors that escape early lead to high delivery. Supranuclear vectors have a much greater chance of success by virtue of their advantageous location. The relative number of supranuclear and cytoplasmic vectors depend on the ratio of nuclear area to the cell area. For cells with this ratio >0.06 (which is true for most cells), the contribution of supranuclear vectors is dominant. Hence, high rate of escape and unpacking maximizes efficiency (Fig. 9). Standard PEI25kDa polyplexes, however, are far away from this maximal delivery state. The calculations show that the ratio of the maximum efficiency 
to the base state (
) is 80. Note that this enhancement ratio is based only on optimization of the intracellular pathway. Other factors such as toxicity of the polymer are not considered here.
Global optimization of vector trafficking in vivo
In contrast to in vitro situations, cells under in vivo conditions are highly three-dimensional (24) and of a stellate shape, with a large central spherical portion (Fig. 1 b). Under these conditions, no region of the cell membrane is particularly close to the nuclear membrane. Due to this major difference in cellular morphology, the spatiotemporal distribution of vectors and hence optimal polyplex parameters are radically different from that in vitro. Fig. 10 depicts the predicted the spatial distribution of polyplexes in three-dimensional cells at 10 min, 4 h, and 24 h post-transfection, calculated using the assumption that the principles and the parameters underlying endocytic trafficking are the same as those in vitro. At short times, there are no vectors within 2 µm of the nucleus. All vectors are located close to the cell membrane, far-away from the nucleus. They are then gradually delivered to the perinuclear region by means of MT-dependent transport of endosomes, and endosome-to-lysosome transfer.
View larger version (20K):
[in this window]
[in a new window]
|
FIGURE 10 Spatiotemporal distribution of vectors in three-dimensional cells. Predicted spatial distribution of PEI25kDa-DNA complexes in a three-dimensional cell at 10 min (solid), 4 h (shaded), and 24 h (dotted) post-transfection.  is the probability that a vector is located at distance r from the nucleus at time t in a three-dimensional cell.
|
|
Fig. 11 shows optimization of vector properties for in vivo gene delivery. Optimal efficiency is observed at intermediate values of transfer and escape rate. The ratio of the maximum efficiency 
to the base state (
) is 40. Early escape of vectors from endosomes leads to underutilization of the endocytic trafficking and consequently, reduces perinuclear accumulation of the vectors, whereas a prolonged delay of vector escape promotes lysosomal degradation of the vector. Similarly, fast unpacking causes premature degradation of DNA inside lysosomes and cytoplasm, whereas slow unpacking limits the amount of DNA available for nuclear translocation. Interestingly, the optimal values of kescape and kunpack are functions of ktr. Thus, to arrive at the optimum pathway, the rates of unpacking and escape cannot be tuned independently but are constrained by ktr.
Impact of polymer degradability on gene delivery
Another important factor in the overall delivery efficiency is the stability of the polymer itself. The foregoing analysis is applicable to nonbiodegradable vectors like PEI25kDa, which are immune to hydrolases, proteases, and nucleases, and offer a great degree of protection to their genetic cargo. For these vectors, it is reasonable to assume that the degradation of DNA is possible only after the vector unpacks. The stability of DNA in both the lysosomes and the cytoplasm is then intimately coupled to the rate of vector unpacking. On the other hand, biodegradable polymers such as poly(ß-amino esters) can themselves undergo hydrolysis in the lysosomes or the cytoplasm. The DNA can thus be degraded without the need to unpack. In other words, there is no direct coupling between unpacking and DNA lysosomal degradation. Surprisingly, this latter case, which makes the DNA more susceptible to enzymatic attacks, leads to an optimal pathway slightly more efficient than that for nonbiodegradable polymers (0.6% vs. 0.4% for in vivo applications). This higher efficiency stems from the inherent flexibility afforded to design of biodegradable polyplexes. Since DNA degradation in lysosomes is no longer conditional upon unpacking of the polyplex, it can be independently optimized as another parameter. Under these circumstances, the efficiency varies linearly with the rate of unpacking kunpack and rate of degradation in lysosomes kdegL, and no optima exist with respect to these parameters (Supplementary Material, section 2.5).
Impact of cell shape on gene delivery
Fig. 12 shows the effects of the shape and the size of in vitro cells on . The shape of the cell is characterized by its circularity and size. The circularity C is defined as C = 4
(A/P2), where A is the area of the cell and P is the perimeter. A small value of C indicates that the cell is elongated. The size is defined as the average distance between the cell membrane and the nuclear membrane. Elongated (less circular) and smaller cells have a larger fraction of the cell area nearer to the nucleus, leading to higher perinuclear concentration and hence a greater delivery efficiency. It is interesting to note that delivery efficiency can vary by almost an order of magnitude, only due to variation of shape and size of the cell.
|
DISCUSSION
|
|---|
Existing models versus current approach
Traditionally, gene delivery by synthetic vectors has been studied and modeled as biochemical reactions between the vectors and cellular components. With recent studies based on single-particle tracking, it has become evident that spatially heterogeneous transport processes are intimately involved in every step along the delivery pathway (21,25). However, existing mathematical models of synthetic gene delivery approximate all transport processes by first-order kinetics (11,12,19). Such simple approximations not only undermine the predictability of the model but also fail to take advantage of data provided by microscopic experiments.
To overcome these issues, we introduced three new features in the present model. First, we introduced spatial coordinates to represent the cell geometry. The cell geometry is reconstructed directly from experimental images. Second, we follow the movements of single particles, just like in single-particle tracking experiments. Thus, intracellular trafficking of polyplexes is no longer represented as discrete jumps from one compartment to another, but as continuous movements of single particles. Third, we use stochastic algorithms to update the positions and the states of polyplexes and cellular organization. This is to reflect the inherent randomness of intracellular processes.
This approach has several advantages. First, it allows us to directly compare model predictions to experimental data. This proves to be crucial in not only model validation but also in verification of our understanding of underlying processes. Second, it enables us to develop a realistic and mechanistic description of intracellular transport phenomena. For instance, it accurately captures the basic physics of microtubule-dependent transport. The particles are allowed to switch intermittently between diffusion and motor-driven directional movements on MTs, just like what they do in experiments. Although the model requires several more parameters to characterize binding and detachment of particles to MTs, these parameters are real, have a physical meaning, and can be estimated directly from microscopic single-particle tracking experiments. Third, the stochastic algorithms grant us a great degree of flexibility in model implementation and solution methods, thus facilitating simultaneous representation of diverse physical and chemical processes that occur at multiple length and timescales. PDE-based models can describe spatially variable processes but are only restricted to simple and fixed cellular configurations. To obtain the solutions to the transport equations for complex geometry and dynamic environment, stochastic simulations are far more superior. Fourth, the model provides an integrated framework to systematically study the influence of all model parameters and extract optimal parameter sets.
The model relates molecular-level binding and trafficking events to whole-cell distribution of polyplexes, and eventually to gene delivery efficiency. The mechanistic treatment of transport allows us to study the effects of cell geometry, motor-assisted transport and endocytic trafficking. For the first time, the spatiotemporal distribution of polyplexes before escape from endocytic vesicles can be predicted with great accuracy. This information is crucial to the design of better vectors.
Perhaps the most important feature of the model is that it represents a step toward a more realistic description of intracellular transport in general and gene delivery in particular. Many processes in gene delivery pathway—for example, endosomal escape and DNA degradation—are extremely difficult to measure experimentally in situ. Accordingly, mechanistically sound models provide a valuable tool in studying the impact of these processes on gene delivery. Such models also provide a logical platform for synthesis and integration of information from diverse sources, bridging scales, cell types, and operating conditions, and allow predictions of cellular and subcellular processes that are inaccessible by current experimental tools.
New insights into gene delivery brought out by this model
The model presented here brings out several novel insights into mechanisms of gene delivery by polyplexes.
The model reveals, quantitatively, that the morphology of cells (e.g., size, shape, and dimension), which has been largely ignored in gene delivery research, strongly influences the delivery efficiency. This is significant since human cells exhibit a wide range of morphologies—muscle cells in the body are very flat and elongated, whereas dermal fibroblasts are stellate and round. It can be inferred from the model results that the strategies to be used for these two applications must be very different. In addition to the morphology of the cell, other cell-specific properties such as endocytic trafficking and microtubule-dependent
TOP
ABSTRACT
INTRODUCTION
