Supporting Text for “Understanding Intracellular Transport ProcessesPertinent 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, CA93106
Contents
Part 1Background Information, Materials and Methods
1.1 SimCell: Simulation of Nanoscale Transport in Cell
1.2 Review of Past Models
1.3 Physical Description of Gene Delivery by PEI-DNA Vectors
1.4 Parameter Estimation
1.5 Reconstruction of Cell Geometry
1.6 Experimental Methods
1.7. Diffusion-advection-reaction equations for the one-dimensional case
Part2Additional Results and Discussion
2.1 Sensitivity Analysis
2.2 Effects of Cell Shape on In Vitro Efficiency
2.3 Comparisons between Dividing and Non-dividing Cells
2.4 Enhancement of Delivery Efficiency by Motor-Assisted Transport
2.5 Comparisons of Gene Transfection by PEI Derivatives to Model Predictions
Part 1: Background Information, Materials and Methods
1.1. SimCell: Simulation of Nanoscale Transport in Cell
The basic modeling objectives of the modeling framework used in the present work are:
- Discrete and mechanistic description of physical and transport processes.
- Realistic representation of cell organelles and components.
- Spatial-temporal evolution of cell as a complex and dynamic system.
We consider two levels of abstraction. The first level of abstraction focuses on the whole-cell level dynamic organization of the cellular components that are involved in the delivery pathways. The second level of abstraction focuses on the interactions of drug carriers with cellular components. Currently, we assume a weak feedback loop between nanoscale carriers and cellular function. This assumption is reasonable since nanoscale carriers generally behave like tracers within the cellular environment, i.e. their presence does not significantly alter or disrupt cellular functions. In future, the effects of perturbations caused by nanoscale carriers to the cellular system can be incorporated into the model.
Mathematically, we consider a state variable , which represents the cellular environment and a state variable , which represents the physical location and biological compartmentalization of the gene vectors.The evolution of is described by a probabilistic nonlinear operator as follows(1, 2),
(1.a)
where includes parameters and constitutive relationships. Similarly, the state of the gene vectors is updated by ,
(1.b)
where consists of parameters that characterize the physicochemical properties of the gene vectors and their interactions with the cellular environment.
The model consists of four main modules (see Figure S1):
(a) Cell-World, which contains cell-specific information, such as cell geometry, distribution of cellular organelles, endocytosis rate, transport properties of endosomes and lysosomes, the number of endosomes and lysosomes, enzymatic activities of endosomes and lysosomes.
Geometry Builder is a standalone application that permits 2D and 3D reconstruction of cellular topology based on experimental images.
(b) Carrier-Physicochemistry, which contains information on carriers and the parameters that govern their interactions with cellular components, such as rates of endosomal escape, unpacking, diffusion coefficient. This module also contains information of the delivered DNA, e.g. half-life of DNA inside the cytoplasm, rate of nuclear entry, etc.
(c) Cell-Physics, which contains mechanistic models and basic numerical algorithms
(d) Cell-Solution, which contains the numerical algorithms that update the states of the system.
In other words, Cell-World and Carrier-Physicochemistry store information on the model parameters and and state variablesand ,Cell-Physicsprovides the numerical realization of the nonlinear operator and and Cell-Solution is the solution of the evolutionary equations A.1.a and A.1.b.
Figure S1. Organization of the computational framework.
1.1.1Mathematical Representation of Cell Geometry: For in vitro applications, we employ a 2D representation of cells. This imitates a cell spread on a surface. Nanoscale carriers bind to and enter the cell from the top surface. The geometry of the cell, including the locations and shapes of subcellular organelles (e.g. nucleus), are reconstructed from microscope images. We create a library of 541 cell geometries based on phase contrast images of ~600 cells using methods described in section 1.6. The cell geometry defines the computational domain. For most calculations, the cell geometry is unchanged during the course of the simulation. This is motivated by the fact that cell motility occurs at a much slower time scale in comparison to other physical processes and therefore has negligible effects on the simulation results. Viewing from the top, the cell interior is divided into two regions, cytoplasmic and supranuclear, separated by the nuclear boundary (see Figure S2). Due to deficiency or complete lack of microtubules in the nuclear region, as seen in various visualizations of microtubule networks, the main transport mode in the supranuclear region is diffusion (3). This is also confirmed by tracking of movements of endosomes and lysosomes in the supranuclear region.
For in vivo applications, we reconstruct the cells from histological and electron microscopy images of fibroblasts embedded in tissues(4). Our reconstructed 3D cells are in good agreement with past descriptions of fibroblasts under in vivo conditions (5, 6).
Figure S2. Schematic of a 2D in vitro cell.
1.1.2.Mathematical Representation of Microtubules: MTs are modeled as straight lines, which grow radially in random directions from the microtubule organizing center(MTOC) towards the cell periphery. MT density is estimated from immunolabeling experiments (described in section 1.4). MTOC is placed randomly in the vicinity of the nucleus. To describe microtubule polymerization dynamics, we use stochastic equations which follow the individual MT length history (7-9). A microtubule switches intermittently between 3 phases, growth, shrinkageand collapse, as follows:
(2)
where q and m are random numbers between 0 and 1. kgrow , kshrinkand kcat are the reaction rate constants for MT growth, shrinkage and collapse, respectively.
1.1.3. Mathematical Representation of Endosomes: The movements of endosomes are approximated as stochastic trajectories of discrete particles of finite sizes. Conceptually, the vesicles are viewed as discrete particles that switch stochastically between free diffusion in cytosol and directed movements on MTs (3, 8, 10-12). To represent these distinguishable physical processes, we introduce substate variable s, s = 1 and 0(0: diffusion, +1: plus-end directional transport and -1: minus-end directional transport). We defineand as velocities of directional movements toward plus-end and minus-end of microtubules. and are first-order rate constants characterizing particle’s transition from s=0 to s= +1 and s=-1, respectively. and are the rate constants of corresponding reverse processes. Thus, the cytoplasmic transport pattern of endosomes is represented by 7 parameters, D0,e(free diffusion coefficient),,, ,,and . The rate constants are lumped representations of several biophysical processes including complex vesicle-motor-microtubule interactions. Since these parameters are not available in literature, we perform independent experiments based on multiple-particle tracking techniques to obtain them (described in section 1.4).
1.1.4. Mathematical Representation of Lysosomes: Lysosomes are also approximated as discrete particles. Unlike endosomes, lysosomes form clusters in the perinuclear region, which are generally immobile. Thus, lysosomes switch stochastically between two different movement modes. When they are free (not inside any clusters), they are capable of recruiting molecular motors and participate in microtubule-dependent transport just like endosomes(13). Therefore, we use the same description for endosomes to describe movements of mobile lysosomes. When lysosomes are inside clusters, they exhibit minimal diffusion. Switching between free and clustering modes is controlled by the rate of clustering, , and the rate of de-clustering . In particular, the rate of clustering is space-dependent and proportional to the local concentration of lysosomes, which is directly measured from single-particle tracking experiments (described in section 1.4). Visualization of formation and breakup of lysosome clusters by fluorescence imaging enables us to quantify transitions of lysosomes between free and immobile states (within a cluster).
Thus, the cytoplasmic pattern of a lysosome is represented by 10 parameters: 7 for free lysosomes, (Dl,f,,, ,,and ), 1 for clustering lysosomes (Dl,c) and 2 for the rates of clustering and de-clustering (and ).
1.1.5.Mathematical Representation of Drug Carriers. The intracellular pathway of a drug carrier can be characterized by five distinct biophysical states : bound to cell membrane (M), inside endosomes (E), inside lysosomes (L), inside cytoplasm (C), and unpacked plasmids (P). Each biophysical state has a distinct transport pattern. If or , the location of the nanoscale carrier is determined by the endosome or the lysosome carrying it. For other states (), transport properties of a gene vector are controlled by its ability to recruit molecular motors, the local viscosity and the presence of obstacles. Transition from one biophysical state to another depends on various biological and physical factors, such as location, interactions with cellular organelles, and molecular components of gene vectors. In the present model, we approximate transitions between the biophysical states of the vectors as first-order reaction, see below. For instance, kescape is the endosomal escape rate, which characterizes the transition from to , and ktr is the transfer rate from endosomes to lysosomes, which characterizes the transition from to .Note that the physical location of the vector remains the same following transitions of biophysical states.
At time t=0, the nanoscale carriers are randomly bound to the cell membrane. They are then delivered to endosomes with an internalization rate, kint. The movements of vector-carrying endosomes are simulated by the transport laws prescribed in section 1.1.3. The probability of transfer to lysosomes is determined by the global transfer rate, ktr, and the local concentration of lysosomes clyso,i. clyso,i is evaluated by counting how many lysosomes are within a radius rl . rl = 3 m provides a robust estimation of clyso,i., which is updated every 20 time steps. The instantaneous, and local endo-to-lyso transfer rate for a vector i is , where is the average lysosome concentration. Following transfer, the movements of vector-carrying lysosomes are updated according to the algorithms given in section 1.1.4. Other transport and reaction events, including diffusion, and unpacking of cytoplasmic vectors, enzymatic degradation and nuclear binding of plasmids are captured using the same approach.
1.1.6. Mathematical Representation of Reaction Events.We approximate transition (reaction) events, such as binding to and unbinding from MTs, and transitions between biophysical states of gene vectors, by Markov processes. A Markov process is a mathematical idealization in which the future state of a system/component is affected by its current state but is independent of its past. This type of memory can be modeled in terms of a transition probability. For example, consider a simple continuous-time Markov process of a random variable S representing the transition between stage A and stage B
(3)
The probability, given S =A at time t, of an AB transition occurring in the interval [t, t+t] is
(4)
We compare a random number uniformly distributed between 0 and 1 to determine whether the event occurs in the time interval t or not.
For instance, endocytosis of membrane-bound vectors is modeled by changing the biophysical state of the nanoscale carrier from membrane-bound to () to caged-inside-endosome (). The probability that such event occurs in the interval [t, t+t] is, where is the rate of internalization.
1.1.7. Mathematical Representation of Transport Events. The position of a diffusive entity (s = 0) is updated using two-dimensional random walk. For MT-bound entities, the following equations of motion are applied,
(5)
where is the position of entity i and characterizes the microtubule j that the entity i is associated with.
Endosomes are capable of alternating between diffusion and MT-dependent transport. If a diffusing endosome is within a radius re from a neighboring MT, the probability that it will bind to the MT is determined by the rate of binding, k. Similarly, detachment from MTs is determined by the rate of unbinding, k’. A similar approach is used to update movements of lysosomes.
1.1.8Simulation Scheme. 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 a library of 500 different cell geometries, which defines the computational domain. Each cell contains 100 endosomes, 100 lysosomes, 1200-1500 microtubules, and 50 gene vectors. We update the states of the system according to the following order: (i) cell geometry, (ii) microtubules, (iii) endosomes, (iv) lysosomes, (v) targets (if necessary) and finally (vi) nanoscale carriers.
To guarantee accuracy and stability of the numerical scheme, we use a small time step, t = 1 s. The time step is much smaller than the governing time scales of the transport processes. Varying t between 0.1s and 2 s does not significantly influence the model predictions.
1.19Data Analysis. Forin vitro simulations, we sample 10000 trajectories from 50 cells. The delivery efficiency of non-dividing cells is calculated based on the number of DNA molecules that gain entry to the nucleus at a final time Tf, at which we terminate the simulations. To obtain the delivery efficiency for dividing cells, we measure the number of plasmids already inside the nucleus and the number of vectors bound to the nuclear membrane, at the time of mitosis.
Alternatively, we can evaluate the delivery efficiency step-by-step. First, we define as the probability that a PEI25kDa-DNA complex, having escaped from an endocytic vesicle at distance re from the nuclear boundary at time , will successfully deliver DNA to the nucleus of a non-dividing cell at a time . The longer the post-escape time , the more likely is the vector to find the nucleus and complete delivery. Second, we introduce, i.e. the probability that a vector, released from endocytic vesicles at time , will successfully complete DNA delivery at some final time . Mathematically,, where is the maximum distance to the nuclear boundary, and is the spatial distribution of vectors inside endocytic vesicles at time te. and can be accurately predicted using simulations. In other words, is a spatially-averaged representation of the local success probability P1. To emulate experimental protocols, Tf is taken to be the time at which cells typically are assayed for gene expression (Tf = 24hrs, typically). Finally, the overall delivery efficiency for in vitro cells can then be calculated as , where is the distribution of escape time. Since, we model escape as a first-order rate process, is an exponential distribution with a mean . Efficiencies are defined in a similar manner for in vivo cells.
1.2. Review of past models.
The kinetic models provided a major boost to the field by providing a mathematical description of intracellular drug delivery for the first time (14-21). The basic assumption of the kinetic models is that all interactions and transport processes can be treated as chemical reactions, characterized by chemical species and kinetic rate constants. The concentration of the reacting species in different cellular compartments as a function of time can be obtained by solving ordinary differential equations (ODEs). Kinetic constants are obtained from qualitative experimental observations, theoretical considerations, and mostly from a direct fit of experimental data. The kinetic models, though simple and effective in representing macroscopic data, do not encapsulate sufficient complexity to enable a meaningful interpretation of all model parameters. The major shortcoming of kinetic models is simplified description of spatial organization and transport processes. For instance, most authors used a simple first-order kinetic equation to describe transport of cytosolic gene vectors to the nucleus: , where and are the number of vectors in the cytosol and in the nucleus (15-17, 19). The nuclear transport rate, kn, is a lumped representation of several processes, including MT-based trafficking, diffusion, and binding to and translocation through the nuclear envelope. The numerical value of kn was obtained by fitting the model predictions to experimental data. This approximation does not reflect the fact that the rate of nuclear entry depends on the local concentration of vectors bound to the nuclear envelope. Accordingly, kn does not have a real physical meaning, rather it is adjusted to match the temporal evolution of or . Fitting of kn to experimental data also involves great uncertainties. This leads to inaccurate predictions of nuclear localization and gene expression outside the range conditions that were used to fit the model.
Other inadequacies arise from the ensemble-average nature of the kinetic approach. The model replaces the individual microscopic identities of drug/gene carriers and cell organelles by macroscopic concepts of “concentration” and “compartment”. Therefore, several important effects can not be modeled, such as spatial distribution, influence of structural and morphological properties of cell organelles, cooperative/competitive effects and fluctuations due to finite number of entities. Furthermore, the interactions and processes involved in intracellular drug delivery spans several scales. Therefore, a macroscopic description of the delivery pathway can not take advantage of information obtained at microscopic and molecular levels. Moreover, macroscopic models can not be used to extrapolate results obtained under in vitro conditions to in vivo conditions. In other words, due to recent advances in experimental techniques, the data available on delivery of nanoscale carriers is much richer than what can be captured by pharmacokinetic models. The modeling challenge, hence, is to sufficiently represent the complex, multi-scale, and dynamic nature of cell. The objective of this study is to formulate a model that captures these inherent characteristics of cellular processes with the emphasis placed on those relevant to intracellular delivery of nanoscale drug/gene carriers. Such holistic approach provides a new paradigm in cell biology (22, 23).
1.3. Physical Description of Gene Delivery by PEI-DNA complexes
Background: Over the past decade, synthetic gene vectors have emerged as a viable strategy for gene therapy(24). Synthetic vectors alleviate safety concerns associated with viral vectors; however, they suffer from extremely low transfection efficiencies. Even for polyethylenimine (PEI)-DNA complexes, which are considered as the gold standard of synthetic vectors, the probability of reaching the targeted cell nucleus is only ~10-4-10-3(25, 26). The low transfection efficiency, in turn, arises from inefficiencies of one or more biophysical processes that determine the intracellular trajectories of the DNA complexes.