Keywords

1 Introduction

Within the past two decades, numerous improvements have been made in the in-vivo and in-situ visualization and characterization of transport processes in cells. Techniques such as confocal microscopy, stimulated emission depletion (STED) microscopy [50, 118], dual color localization microscopy (2CLM) [42], spectral position determination microscopy (SPDM) [76] and enhanced staining [45, 108] allow for detailed and statistical analysis of translocation trajectories and distribution of intracellular organelles or marked enzymes and enzyme complexes down to the nanometer scale [26]. The number of spatiotemporally resolvable and quantifiable intracellular transport processes is thus constantly increasing and the complexity of the resulting network of cooperating and competing transportation events is rising. Interaction between different cells such as clustering and quorum sensing adds a further dimension of complexity.

The sheer complexity of such intra- and intercellular biological processes increasingly drives the need for the development of sophisticated modeling techniques in order to understand causalities and provide tools for the determination of bottlenecks and optimization of biotechnological processes. The increasing amount of available data combined with still exponentially growing computation power allows speculation about the possibility of whole-cell simulations, i.e. the construction of virtual cells in silico, in the not-so-remote future, an idea that has been especially in vogue at the beginning of the millenium [130]. Computational implications of whole-cell simulation with respect to the complexity of intracellular proteome and metabolome diversity, as analyzed in [119] and [6], suggest that it may in principle be possible to perform a cell cycle of, e.g., E. coli within a couple of weeks on powerful computer clusters, provided that nanoscale (i.e. molecular) effects are widely neglected. However, this holds only true if parallelization and therefore high scalability can be provided, which is a non-trivial problem [6, 92].

In recent years, a multitude of specialized models for well-defined biological questions have been published. Meanwhile, more generalized and multi-purpose software packages have been developed in order to ease handling for non-specialists. This chapter will address both aspects, with the focus on transportation processes within cells. We start with simple free-diffusion compartmented models in Sect. 2, then proceed to active intracellular transport (Sect. 3), to quantitative descriptions of vesicle distribution (Sect. 4), and membrane transport (Sect. 5). We consider the influence of cell morphology fluctuations and cell growth, as well as the impact of transport variations during the cell cycle in unsynchronized cultures (Sect. 6), and the dynamics of protein translation (Sect. 7). In the following section, the issue of standardized model data interchange will be addressed, and finally various general-use software packages for cell and/or biochemical reaction network simulation will be introduced.

2 Free Diffusion

Biochemical models, in their simplest form, assume that every compartment under consideration is perfectly mixed and hence neglect transport processes exhibited by the contributors to the biochemical reactions modeled . This can be a reasonable approximation if transport and diffusion processes are considered to be fast compared to the actual reactions. In many cases, however, due to the relatively large volume of eukaryotic cells and their complex compartmentation, the spatiotemporal motion dynamics of substrates, reaction products, co-factors, vesicles or even the reaction sites themselves (e.g. fluctuating membrane-bound protein complexes) cannot be neglected. For this purpose, it is necessary to incorporate a reasonably resolved three-dimensional view of the cells and/or compartments under investigation, e.g. obtained using confocal electron or light microscopy. The numerical simulation of diffusion processes (Wiener process, Fick’s law) in predefined three-dimensional geometries is within the scope of standard cell simulation packages, which will be summarized in Sect. 9.

The limiting factor in free diffusion processes is denoted by the diffusion coefficient in the cytosol, \(D_{\rm cyto},\)which can be given either in absolute units (e.g. \(\upmu \hbox{m}^2/\hbox{s}\)) or as relative diffusion coefficient compared to water \(D_{\rm cyto}/D_{\rm water}.\) The diffusion coefficient strongly depends on the properties of the solute and solvent. In general, this is defined by the Navier–Stokes equation:

$$ D=\frac{k_B T}{6 \pi \eta R_0}, $$
(1)

with \(k_B\) denoting the Boltzmann constant (\(k_B = 1.380 6504\times 10^{-23} \hbox{J/K}\)[93]), T the temperature (in K), \(\eta\) the dynamic viscosity of the solvent, and \(R_0\) the hydrodynamic radius, or Stokes radius. In real processes, the Stokes radius is roughly correlated with the molar mass [10], leading to generally rapid diffusion of small solutes compared to solutes of high molecular weight. Due to their high diffusion rates, small molecules or ions are therefore suitable to act as messengers in intracellular signalling networks, which have been reviewed in [71].

The diffusion coefficients for numerous substances have been determined empirically, usually by application of fluorescence recovery after photobleaching (FRAP) [112, 129]. Experimentally derived diffusion coefficients of substrates of different type and molar weight in numerous compartments have been reviewed in [126]. It is noted that small solutes (\(\hbox{BCECF, MW} = 550-820\, \hbox{Da}\) [61]) diffuse approximately four times slower in the cytosol than in water, a finding similar to that of FITC-dextrans and FITC-Ficolls up to a \(\hbox{MW}\approx500\,\hbox{kDa},\)green fluorescent protein (GFP)\(( \hbox{MW}\approx27\,\hbox{kDa},\,D_{\rm GFP, cyto}=2.5-3\times 10^{-7}\hbox{cm}^2 s^{-1})\) and very small DNA fragments of not more than\(\approx 100\,\hbox{bp},\)i.e. \(\hbox{MW}\approx70 \,\hbox{kDa}.\) In any case, whenever the limiting molar mass of a specific substrate is exceeded, the diffusion rate decreases in a highly non-linear fashion to almost zero [83], an effect which is presumably mainly determined by the influence of the actin cytoskeleton [19]. Diffusion measurements of unconjugated GFP in mitochondria \((D_{\rm GFP, mito}=2-3\times 10^{-7}\,\hbox{cm}^2\,\hbox{s}^{-1})\) [106] indicate similar free diffusion characteristics, while the apparent diffusion in the endoplasmic reticulum (ER) is decreased approximately three-fold \((D_{\rm GFP, ER}=0.5-1\times 10^{-7}\,\hbox{cm}^2 \,\hbox{s}^{-1})\)[21], which is suspected to be caused mainly by the complex geometry of ER.

The a-priori prediction of free protein diffusion kinetics in aqueous solutions according to their structure and conformation has been investigated in the past. A numerical prediction model for protein diffusion coefficients in water was proposed in [11]. It can be adapted to a wide range of protein surfaces. The model was compared to conventional approaches of Einstein and [36] and yielded comparably superior results with respect to experimental data of lysozyme and tobacco mosaic virus. A more recent model based on the Yukawa potential was published in [43].

The estimation of protein diffusion on charged membranes was the focus of the work published in [70]. This model considers lipid lateral reorganization and demixing when adsorbed charged macromolecules diffuse on membranes in a time scale of microseconds, which is several orders of magnitudes faster than conventional molecular dynamics (MD) approaches.

3 Active Transport

In eukaryotic cells, a significant portion of the intracellular transport of particles and vesicles is facilitated by active transport along microtubules (MT) and/or actin filaments (AF). This is necessary due to the size of eukaryotic cells and their complicated compartmentation, which would make particle transport by pure diffusion very inefficient. The resulting process has similarities with “facilitated diffusion” which can be found in other biological processes such as oxygen transport by myoglobin [142, 143] and hemoglobin [88], rapid association of DNA binding proteins with specific DNA sequences [47, 66], or ribosomal subunit translocation [65, 69]. MT are polymers that are mostly formed by nucleation starting on microtubule organization centers (MTOC) [82]. They exhibit a characteristic fluctuation pattern throughout the cell cycle [100], usually distinguished by growth, shrinking, and pause phases; and additionally occasional “catastrophe” events [16].

Simple compartmented models as described in Sect. 2 assume force-free diffusion of particles and substances within each compartment. Hence, directed transport is neglected. In recent years, increasing effort has been put into the mathematical description of intracellular active transport processes. Special attention has been laid on the non-continuous, but saltatory active movement of particles (“random walk”), i.e. the fact that particles are usually transported for distances of up to \(5-20 \, \upmu\,\hbox{m}\) at more or less constant velocities, with intermediate pauses of more than 1 s and even occasional reverse transport, depending on the region within the cell [72]. An early publication [128] assumed different transport velocities (fast anterograde, slow anterograde and retrograde), which were later questioned [3]. The work of [123], based on improved intracellular movement classification approaches [139], describes an analytical mathematical and one-dimensional model of motor-assisted transport via MTs or AFs of intracellular particles. It distinguishes between an unidirectional model with all filaments having the same polarity, and a bi-directional model. It incorporates simple kinetics for the association and dissociation of particles to and from MTs and/or AFs. It allows for free diffusion of unattached particles and constant motion of attached particles based on partial differential equations. The motion of particles satisfies the one-dimensional transport equation

$$ \begin{aligned} \frac{\partial n_0\left(x,t\right)}{\partial t} -D\frac{\partial^2n_0\left(x,t\right)}{\partial x^2} &= -\left(k_++k_-\right)n_0+k_+^{\prime}n_++k_-^{\prime}n_-\\ \frac{\partial n_\pm\left(x,t\right)}{\partial t} + v_\pm\frac{\partial n_\pm\left(x,t\right)}{\partial x} &= k_\pm n_0 - k_\pm^{\prime}n_\pm, \end{aligned} $$
(2)

with \(n_0\left(x,t\right)\) denoting the density of free particles at distance x along the filament at time t, and \(N_\pm\,\left(x,t\right)\) the densities on right- and left-directed filaments. \(k_{+}\) and \(k_{-}\) are the first-order rate constants for binding to filaments, and \(k_{+}^{\prime}\) and \(k_{-}^{\prime}\) for detachment, respectively. The motor velocities are denoted by \(v_{+}\) and \(v_{-}.\) The particle flux \(J\left(x,t\right)\) consequently reduces to

$$ J\left(x,t\right)=-D\frac{\partial n_0\left(x,t\right)}{\partial x}+n_+\left(x,t\right)v_++n_-\left(x,t\right)v_-. $$
(3)

From these equations, one can derive the mean free-diffusion lifetime \(\tau_{\rm off}=\left(k_++k_-\right)^{-1},\)the mean free-diffusion path length \(l_{\rm off}=\sqrt{D/\left(k_++k_-\right)},\)the mean active transport lifetime \(\tau_\pm=1/k_\pm^{\prime},\)and the mean active transport path length \(l_\pm=\left|v_\pm\right|/k_\pm^{\prime}.\) Under the assumption that free diffusion is absent, the mean drift velocity due to active transport is equal to

$$ \overline{v}=\frac{K_+v_++K_-v_-}{K_++K_-+1}, $$
(4)

with \(K_\pm\equiv k_\pm/k_\pm^{\prime}.\) The effective diffusion constant yields

$$ D_{\ast}=\frac{D+K_+\left(v_+-\overline v\right)^2/k_{+}^{\prime}+K_-\left(v_--\overline v\right)^2/k_{-}^{\prime}}{K_++K_-+1}. $$
(5)

The model of [123] considered diffusion as a one-dimensional process along each MT/AF. If a rotationally symmetric distribution of MTs/AFs in the cell is assumed, the extension to three dimensions is relatively straightforward. In the work of [34], the model has been extended with respect to three-dimensional diffusion in a cylindrical neighborhood of each MT or AF.

This, however, seems to show that pure analytical mathematical models are not appropriate to cope with complex (i.e. non-symmetric and non-periodic) cell, MT, and AF geometries, since the number of equations to solve would be hard to handle when dealing with variable MT and AT construction and destruction kinetics, non-linear MT and AF geometries, arbitrarily-shaped compartments etc. Thus, numerical and hence spatiotemporally quantized models based on Monte Carlo methods are becoming more prevalent in order to obtain more flexibility with respect to the aforementioned model properties.

A computational model for quantification and analysis of the switching process between MT and AF transport has been published by [122]. This model simulates pigment transport in melanophores. Here, radial particle motion is classified as three discrete states—movement to the MT plus end \((P_{+}),\) to the MT minus end \((P_{-}),\) and pause\((P_0),\) with the corresponding transition rate constants \(k_{1-6}.\) Transport on AFs is considered as two-dimensional diffusion, and the transition rates for the switching of a particle from MTs to AFs or AFs to MTs are denoted by \(k_{\rm MA}\) and \(k_{\rm AM},\) respectively. The cell geometry is idealized to a 2D circular shape and MTs are considered to form an ideal radial array. The resulting set of linear differential equations has been solved numerically, since analytical treatment turned out to be impractical due to the fact that the unknown parameters enter both the equations and boundary conditions. The numerical solution has been obtained using the VirtualCell framework [121] (see also the introduction of generalized cell simulation toolboxes in Sect. 9) on a mesh size of 0.5–0.55 \(\upmu\hbox{m}.\) The parameter fits found in this study (\(k_{\rm AM}\ll k_{\rm MA}\) during dispersion) suggest that the transfer from MTs to AFs is almost irreversible during dispersion.

Intracellular transport processes for synthetic, i.e. non-viral, gene delivery during transient gene transfection processes have been the subject of a simulation approach published in [23], which is based on preceding works [24, 103, 123]. It utilizes stochastic modeling in order to obtain the highest flexibility with respect to arbitrary cell and nuclear morphology and intracellular organization. No model fitting has been incorporated.

A non-linear and saltatory MT dynamics model is included that utilizes frequent growth and shrinkage cycles, as well as rare catastrophe, i.e. complete MT destruction, events. In this model, the transfected plasmids, mediated by polyplexes such as polyethylenimine (PEI), are internalized via endocytosis. The endosomes generated can be subject to either diffusion or active transport along MTs in the positive or negative direction. The contained plasmids can unpack from the polyplex, underlie degradation, escape to the cytosol or transfer to lysosomes. After leaving vesicles, the plasmids and/or polyplexes can only move passively by diffusion. Cell division and mitosis are neglected, as well as cell growth. The simulation results are compared to measured data of human skin normal fibroblasts.

It has been proven by simulation that the results of in-vitro transfection optimization cannot easily be transferred to in-vivo conditions. For example, the transfection efficiency is highly dependent on the endosomal escape rate of plasmid vectors and/or polyplexes; overshooting the maximum can lead to an up to 1,000-fold reduction in delivery efficiency. It can be concluded that since cell morphology plays a significant role in the transfection process, the stochastic simulation scheme presented, though containing several simplifications, provides a flexible and consistent quantitative description of intracellular transport processes during gene delivery.

More detailed analyses of MT dynamics and transport have been performed in [2], finding that the growth and shrinking behavior of freshly assembled MT differs significantly from those close to the cell membrane. Based on this data, we derived a detailed MT dynamics model that has been implemented in a stochastic model developed in our group for the quantitative analysis of transient gene transfection processes in mammalian cells. The MT model is two-fold: a strongly growth-favoring variant for nascent MTs (phase I) with a length of less than 85% of the cell diameter is implemented as a Gaussian distribution fitted to measured growth velocites (Fig. 1a, with the mean velocity \(v^I_{\upmu}=16.6\,{\upmu\hbox{m}\,\hbox{min}^{-1}}\) and standard deviation \(v^I_{\sigma}=13.0\,{\upmu\hbox{m}\,\hbox{min}^{-1}}.\) For MTs in phase II, i.e. those which have exceeded 85% of the cell diameter, a different distribution can be derived, as shown in Fig. 1b. The number of shrinking MTs is significantly higher, thus two combined Gaussian fits are employed with \(v^{II}_{\alpha,\upmu}=-24.9\,{\upmu\hbox{m}\,\hbox{min}^{-1}}, \,v^{II}_{\alpha,\sigma}=11.4\,{\upmu\hbox{m}\,\hbox{min}^{-1}}\)(shrinking), and \(v^{II}_{\beta,\upmu}=17.4{\upmu\hbox{m}\,\hbox{min}^{-1}}, \,v^{II}_{\beta,\sigma}=5.4{\upmu\hbox{m}\,\hbox{min}^{-1}}.\) The fraction of MTs in rest, which has been measured independently, is denoted by \(k^I_0\) for phase I and \(k^{II}_0\) for phase II. In rare cases, MTs can collapse completely (catastrophe), an event defined by the mean collapse frequency \(f_{\rm c}.\)

Fig. 1
figure 1

Adaptation of microtubule (MT) dynamics for a two-stage model and Gaussian fits. a For nascent MTs. b For MTs near the cell membrane (data adapted from [72])

Full size image

In a simulation framework describing transient transfection dynamics in mammalian cells [63], the impact of active microtubule-facilitated polyplex transport compared to pure passive transport can be simulated depending on the cell geometry (nucleus-to-cell diameter ratio \(\Phi_{\rm nucl}\) and eccentricity of the nucleus \(\epsilon_{\rm nucl}\)), and on clustering of the cells, as visualized in Fig. 2. The transfection efficiency increases by approx 9–10% in the range of \(\Phi_{\rm nucl}=0.7\) and centered nucleus \((\epsilon_{\rm nucl}=0),\) no matter whether cells are clustered or not. With the nucleus being more off-center \((\epsilon_{\rm nucl}\rightarrow 1),\) the improvement decreases to almost \(0\%\) in the case of non-clustered cells, but increases to \(37\%\) in the case of clustered cells. This effect is stronger for smaller nuclei sizes.

Fig. 2
figure 2

Simulated improvement of lipoplex transfection efficiency in HEK293s cells by active transport via MT, depending on nucleus-to-cell size ratio, off-centering of the nucleus and clustering of cells. a Relative increase of with active transport based on MT (without cell clustering). b With cell clustering. c Bottom plane: without clustering; top plane: with cell clustering

Full size image

Further complications when considering active MT transport arise from the fact that MTs usually do not really grow linearly and radially, as usually assumed. A detailed analysis of MT bending characteristics within the microtubule array in 3T3 fibroblasts has been performed in [2], yielding the result that, although the bending angle does usually not exceed \(10-15\,^{\circ}\) within \(5\,{\upmu\hbox{m}},\) the amount of strongly bent MTs rises with distance from the MTOC (60% at \(10\,{\upmu\hbox{m}}\) distance from the MTOC, and 75% at \(20\,{\upmu\hbox{m}}.\)) Another issue to consider is the de-novo formation of non-centrosomal MTs, which are evident in neurons [136], but have also been discovered, e.g., in A498 cells [144].

4 Dynamics of Vesicle Distribution

The spatiotemporal distribution of vesicles, especially transport vesicles, in eukaryotic cells is subject to a complex regulation system that is still not completely understood. A generalized theory of spatial patterns of intracellular organelles has been reported in [22]. It describes the flow of organelles within cells, which is mediated by three transport processes, driven by kinesin, dynein, and myosin. Interactions between organelles are not considered.

The transport processes are denoted by s, with \(s=0\) for free diffusion, \(s=-1\) and \(s=+1\) for microtubule (MT)-mediated transport in the minus or plus direction, respectively, and \(s=2\) for actin filament (AF) transport. The corresponding attachment and detachment rates are \(k_s\) and \(k^{\prime}_s,\)respectively, with the affinity constant being \(K_s=k_s/k^{\prime}_s.\) The organelle density at time t at radial position r is denoted by \(\tilde c_s\left(r,t\right).\)After nondimensionalization and approximation over all states \(\tilde c \left(r,t\right)=\sum{\tilde c_s\left(r,t\right)},\)the model yields

$$ \Pi\frac{\partial c}{d \tau}\approx \Phi\frac{1\partial}{\xi}\frac{\left(c\xi\right)}{\partial\xi} + \Omega\frac{1\partial^2}{\xi}\frac{\left(c\xi\right)}{\partial\xi^2}+ \Delta\frac{1\partial}{\xi\partial\xi}\left[\xi\frac{\partial c}{\partial\xi}\right], $$
(6)

with

$$ \begin{aligned} \xi &=r/R_c, \xi_N=R_N/R_C, \tau=tV/R_C, c=\tilde c/C_0, \\ \Pi &=1+K_{+1}+K_{-1}+K_2, \Phi=K_{-1}-K_{+1}, \Delta=\tilde D_2 K_2+\tilde D_0, \\ \Omega &=\frac{V}{R_C}\left( \frac{K_{+1}}{k^{\prime}_{+1}}\left(1+\frac{\Phi}{\Pi}\right)^2 + \frac{K_{-1}}{k^{\prime}_{-1}}\left(1-\frac{\Phi}{\Pi}\right)^2 + \frac{K_2}{k^{\prime}_2}\left(\frac{\Phi}{\Pi}\right)^2 \right), \\ \tilde D_2 &= D_2 / VR_C {, \hbox{and} }\, \tilde D_0/VR_c. \end{aligned}$$

Directed motions on MTs are represented by \(\Phi>0\) (toward nucleus) and \(\Phi<0\) (away from nucleus), dispersive radial motions along MTs are represented by \(\Omega,\) and dispersive motions along the cell surface as a combination of cytosolic diffusion and AF movement are denoted by \(\Delta.\) The authors define two Péclet numbers to quantify the relative contributions of each type of motion, leading to a certain equilibrium spatial distribution. First, the one-dimensional Péclet number \(P_{e,1{\rm D}}=\Phi/\Omega\) compares directed and diffusive motions on MTs. Secondly, the two-dimensional Péclet number \(P_{e,2{\rm D}}=\Phi/\Delta\) describes the the ratio of directed and radial motion along MTs to diffusive motion along the cell surface. After solving Eq. 6, four distinct limiting patterns can be determined: (1) & aggregation near the cell center; (2) hyperdispersion which denotes the concentation of organelles near the cell periphery; (3) areal dispersion, i.e. uniform distribution over the cell surface area; and (4) radial dispersion, i.e. uniform radial distribution. A regime map can be derived that provides a correlation between the two relevant transport variables, \(P_{e,\rm{2D}}\) and \(P_{e,\rm{2D}}/P_{e,\rm{1D}}=\Omega/\Delta,\) and the resulting equilibrium spatial distribution. The computed patterns have been confirmed by experimental data [22].

5 Membrane Transport

Biomembranes are lipid bilayers with embedded proteins that serve as selectively passable barriers between either the cell’s interior and its exterior or between different compartments within the cell. The transport of solutes through biomembranes is facilitated by one or more of three distinct processes: passive diffusion, facilitated diffusion and/or active transport. The composition of membranes and embedded proteins is however highly specific and dynamic, a necessary prerequisite for the correct trafficking of nutrients, lipids etc. in and out of the cell and into and from the corresponding organelles. For an overview of the functional diversity of membrane transport, the reader is referred to excellent reviews [48, 89].

5.1 Passive Diffusion

Pure lipid bilayers, i.e. synthetic membranes without proteins, enable small molecules to diffuse passively through the membrane. This spontaneous and concentration-balancing process is strongly dependent on molecule size and polarization. Ions are virtually incapable of diffusion through membranes. The passive diffusion velocity is proportional to the concentration deviation, hence the diffusion can be described by a solute-dependent permeability coefficient P [cm/s].

The diffusion rate can therefore be expressed as

$$ \frac{\partial n}{\partial t} = P \, A \, \frac{\partial C}{\partial x}, $$
(7)

with \(\partial n / \partial t\) denoting the diffusion rate in \( \hbox{mol} \,s^{-1},\)A the considered membrane surface area, and \(\partial C / \partial x\) the concentration gradient.

Typical permeability coefficients P vary in the range of, e.g., \(4.0\times10^{-6}{\hbox{cm}\,\hbox{s}^{-1}}\) for urea, \(3.4\times10^{-3}{\hbox{cm}\,\hbox{s}^{-1}}\) for \(\hbox{H}_2\hbox{O},\) and \(2.9\,{\hbox{cm}\,\hbox{s}^{-1}}\) for hydrochloric acid [137].

5.2 Channels and Facilitated Diffusion

In free diffusion processes, the diffusion of small lipophobic or hydrophilic molecules and ions through the membrane is enabled by membrane-bound channel or port proteins with a hydrophilic interior, forming transmembrane pores. This can either be performed without binding of the channel protein to the substrate, resulting in high transport capacities of \(\gg 10^6\) molecules or ions per second per channel (free diffusion) or with binding and substrate-triggered conformation change of the porter protein. Since the transport is passive, its direction is always from high to low substrate concentration, resulting in a decrease of concentration gradient.

Unlike in passive diffusion as expressed in Eq. 7, the transport rate does not increase linearly with increasing concentration gradient, but reaches a saturation level determined by \(v_{\rm max}\) in \(\rm mol/\left(cm^3\, s\right):\)

$$ \frac{\partial n}{\partial t} = \frac{v_{\rm max}}{ 1+K/\frac{\partial C}{\partial x}}, $$
(8)

with K (in \(\hbox{mol}\,\hbox{cm}^{-4}\)) determining the speed of saturation, i.e. the time elapsed until \(1/2 \, v_{\rm max}\) is reached.

Crucial characteristics of channels and ports are their substrate selectivity and their regulation properties, i.e. how the flux is maintained or interrupted depending on external signals. These signals can be mediated by signal molecules, especially ions, phosphorylation/dephosphorylation, voltage variations etc.

5.2.1 Aquaporins

AQP are special channel proteins dedicated to the active transport of water through membranes [12]. They reach a transport capacity of up to \(3\times10^9\) water molecules per seconds per channel [1] while preventing protons from moving through the channel along the water molecule network using the Grotthuss mechanism [86], which is essential for the conservation of electrochemical gradients. The single channel water permeability of AQP1 is \(\approx 4.6-11 \times10^{-14} \hbox{cm}^3/\hbox{s}\) [114, 145].

The related sub-family of aquaglyceroporins are responsible for transport of other small molecules, such as glycerol, ammonia and urea. Furthermore, aquaporins seem to play a pivotal role in cell migration, e.g. of tumor cells [104, 135].

5.3 Active Transport

During active transport processes, the substrate molecules are transported against the concentration gradient, thus consuming energy that is usually derived from the cell’s metabolism. The transport is enabled by carrier proteins. The transportation process can be suppressed via either competitive inhibitors that bind to the transportation site of the carrier enzyme instead of the substrate, or via allosteric inhibitors that bind to elsewhere than the active site of the carrier, forcing it to change its conformation.

If the transport of single substrate molecules is performed without the need for any second substrate (primary transport), the carrier is described as a uniporter. In this case, due to the need for active binding of the transporter protein or cofactors to the substrate and their limited number, the transport velocity characteristics is analogous to that of facilitated transport as denoted in Eq. 9. In the case of coupled transport in coordination with a second substrate, the carrier is referred to as a co-transporter. The coupled transport can either take place in the same direction (symporter) or in the opposite direction (antiporter). In either case, the passive transport of one substrate along its concentration gradient increases entropy and therefore yields energy, which is in consequence used for the active transport of the second substrate against the concentration gradient (secondary active transport; see also [48]). For coupled transports, the flux can be approximated by the sigmoidal function

$$ \frac{\partial n}{\partial t} = \frac{v_{\rm max}}{ 1+\left(K/\frac{\partial C}{\partial x}\right)^h}, $$
(9)

with \(h>1\) denoting the degree of (positive) cooperativity (see also [73]).

5.4 Microdomains Within Membranes

In recent years, the view of the distribution of functional proteins such as transporters along the membrane (especially cell membrane and ER) surface has changed. Formerly regarded as a homogeneous fluid, it is now believed that there are well organized sections of approx 50 nm lateral size with specialized functions, so-called functional rafts [120]. They are formed in a highly dynamic manner with contributions from sphingolipid and cholesterol [46, 78, 97], presumably in coordination with the underlying actin-based cytoskeletal network [111]. Functional rafts are believed to play a pivotal role in protein sorting in the endoplasmic reticulum (ER) and Golgi apparatus [51].

A common model assumption in the morphology of rafts is that they form self-assembling dynamic lipid shells [4, 87], able to bind specific proteins. The characterization of the motility of membrane-bound proteins has moved increasingly into the focus of investigation, featuring confocal and two-photon microscopy as well as single particle tracking (SPT) with sophisticated tracking algorithms (reviewed in [62, 77]). The “hop-diffusion” motio n characteristics of transmembrane proteins when shuffled between microdomains has been reviewed in [20] and [111]. A common but not yet fully conclusive finding is that the proteins are seemingly able to diffuse freely within a small area and sporadically exhibit a “hop” to another region. A corresponding simulation approach therefore implements a two-stage transmembrane protein diffusion model that involves reduced long-range diffusion by introducing barriers while short-range diffusion remains unchanged [134].

6 Cell Cycle and Cell Morphology Dynamics

Most current intracellular transport

modeling approaches neglect the effect of transport variations during the progress of the cell cycle and the impact of cell growth or cell morphology changes during the cell cycle. Since it is clear that cell and compartment morphology as well as intracellular organization have a strong influence on the individual transport processes and their balance [23], it is worthwhile increasing model accuracy by including cell cycle-induced effects.

The dynamic implications of complex cell cycle regulation have been evaluated based on general physiological properties of the most important regulatory modules of the CdK network regulating the cell division cycle, which have been described in a generic computational model [18]. This model aims at a generic representation of yeast, frog cell and human cell cycle regulation. It contains numerous different control loops, which have been tentatively adapted to yeast and mammalian cells by simplification of sub-networks and clustering. By means of bifurcation analysis, the typical oscillatory behavior of the regulating enzymes can be detected and analyzed. The increase in cell mass (and therefore cell volume) during the cell cycle is presumably connected to the cell cycle control [25],Footnote 1which is incorporated in this model by defining the cytosolic synthesis rate of the regulatory proteins A, B, and E as proportional to the cell mass. The mass increase during the cell cycle is assumed to be proportional to the cell size, leading to exponential growth behavior, i.e. \(m\left(t\right)=m_0\,\mathrm e^{t/\upmu}\) with \(0\leq t < \upmu\,\ln\left(2\right).\) This implies an almost linear diameter growth \(d\left(t\right)=d_0\,\mathrm e^{t/\left(3\upmu\right)} \approx d_0\,(1+0.3654 t),\)assuming that growth is isotropic. However, growth is usually non-constant during the progress of the cell cycle (reviewed in [40]), but can at least be considered linear and isotropic in G1 phase.

A specialized and empirical model for the advance of the cell cycle in myeloma cells has been proposed in [79]. If the cell cycle is denoted by \(^c\overline{t}_\phi = (t-^ct_\delta) / ^c\tau,\) with t representing the simulated time point in s, \(^ct_\delta\) the time point of the last division of cell c, and \(^c\tau\) the cell’s cell cycle duration, the cell cycle phase \(\Omega\) with increasing \(^c\overline{t}_\phi\) according to [79] is defined as

$$ \Omega = \left\{\begin{array}{crl}G1; & 0 \leq^{c}\overline{t}_{\phi}<0.50 \\ S; & 0.50 \leq^{c}\overline{t}_{\phi}<0.82\\ G2M; & 0.82\leq^{c}\overline{t}_{\phi}<1 \end{array}\right. .$$
(10)

Using the aforementioned dynamics of the cell cycle and cell growth, the impact and spatiotemporal transport processes can be estimated. In the transient transfection simulation model described in [63], the retention times of transfected plasmids in different compartments of HEK293s cells can be estimated depending on cultivation and transport parameters, with and without growth (but with same average volume). In this model, the volumetric growth is assumed to be exponential during the G1 phase, and linear in S/G2M phases. Figure 3denotes the corresponding trajectories in the different compartments (cytosol and vesicles combined; nucleus). The (peak) number of plasmids delivered to the nucleus reduces by approx 48% when no growth is assumed.

Fig. 3
figure 3

Influence of growth on transfection dynamics in HEK293s cells. Continuous lines with volumetrically linear growth in G1 phases, dotted lines without growth and fixed averaged size. a Percentage of plasmids, both complexed and pure, in cells outside of nucleus (accumulated for vesicles and cytosol). b Percentage of plasmids in nucleus

Full size image

6.1 Cell-Cycle-Dependent Transport Variations

It has been well known for a long time [113] that significant intracellular transport rate variations can occur during the progress of the cell cycle. Many of them may be related to actin polarization as cells enter the S phase. For example, plasmid uptake rates via endocytosis in synchronized D407 cells can vary by up to 60% [85]. In addition, the activity of nuclear import facilitated by the nuclear pore complex (NPC) is altered drastically [84, 132]. Subsequent transport rate variations of macromolecules between cytoplasm and nucleus have been reviewed in [15, 138].

7 Translation and Protein Transport

The protein translocation and secretion process in eukaryotic cells is highly complex, selective and non-linear [74, 107] and is thus a worthwhile subject for detailed quantitative analysis, especially for optimization of protein expression in pharmaceutical protein production cell lines. However, the intracellular post-translational transport of proteins in eukaryotes has not yet been the subject of detailed mechanistic and spatiotemporally resolved models, unlike in bacteria, for which protein traffic is quantitatively better understood [29]. It is thus desirable to obtain a deeper mechanistic and quantitative understanding of translation and & post-translational processes in order shed light on biological questions such as quantification of the impact of mRNA degradation, protein folding and selection, post-translational modification, glycosylation and the secretory pathway on the expression level.

7.1 Overview and Limiting Factors

The protein transport and targeting during and after translation consists of multiple steps which will be briefly mentioned in this section, with each of them potentially influencing the overall protein production ability. For more detailed insight into the numerous processes contributing to protein translocation, targeting and secretion, the reader is referred to excellent reviews [44, 52, 64, 67, 109].

After transcription and mRNA export from the nucleus, the majority of proteins are synthesized in the cytosol with the help of ribosomes which are either suspended in the cytosol or bound to the rough endoplasmic reticulum (ER) or to the nuclear envelope. The translocation into the ER is conducted in either co-translational (mammalian cells) or post-translational fashion (mostly in yeast). The protein remains in an unfolded or loosely folded state and is transferred across the ER membrane through the translocon either completely (soluble proteins with usually amino-terminal, cleavable signal sequences) or partially (with the hydrophilic ends either crossing the membrane or remaining in the cytosol) [109]. The ER-bound proteins are translocated by a passive channel that associates with three different types of partners to drive the force for translocation: with ribosomes (mostly used for secretory proteins and most membrane proteins) for co-translational translocation, with the Sec62/Sec63 complex (a tetrameric membrane-protein complex) or with the luminal chaperone binding immunoglobulin protein (BiP) for post-translational translocation [44, 109, 110]. While in the translocon, the signal peptidases are cleaved off [96], the nascent chain is glycosylated by oligosaccharyltransferase (OST) and co-translational folding is conducted. Assembly of multicomponent complexes may also occur co-translationally in some cases [64]. After release from the ribosome/translocon complex and having almost reached the final conformation, oligomers are assembled in three stages under the presence of three folding enzymes (BOX 2) [30]. Misfolded proteins are identified in a three-stage process featuring molecular folding sensors and chaperones [30] and possibly retrotranslocated into the cytosol via the translocon for degradation in the cytosol by the cytoplasmic proteasome [74, 133].

7.2 FoldEx: ER Protein Export Model

A non-compartmented protein trafficking model called “FoldEx” describing the folding-end export of proteins from the ER has been published in [141]. It provides a quantitative description of the dynamics of the competing pathways that are reponsible for ER import of the unfolded proteins facilitated by the translocon; high-affinity (ATP-bound) and low-affinity binding to chaperones; recognition of correct, incomplete or uncorrect folding; re-folding procedures; and association to export or retranslocation pathways. An overview of the model is given in Fig. 4 and Table 1. The pathways are implemented utilizing Michaelis–Menten kinetics.

Fig. 4
figure 4

Overview of protein export model according to [141]. See text and Table 1 for explanations and details

Full size image
Table 1 Desciption of pathways and rate constants in the protein export model according to [141], as outlined in Fig. 4
Full size table

The authors conclude that the export efficiency of fast-folding (small) proteins mainly depends on the relative activities of the export pathways and ERAD (ER-associated degradation pathway), i.e. the corresponding concentrations ratio \(c_E/c_R,\) while being largely constant with varying \(k_8\) at wide ranges. On the other hand, the export efficiency of slow-folding proteins is more sensitive to folding kinetics; thus even stable proteins may be degraded to a considerable extent when they fold slowly (i.e. transmembrane proteins). Misfolded proteins may re-enter the ER-assisted folding (ERAF) pathway to some extent, increasing the export efficiency of proteins prone to misfolding.

7.3 Protein Targeting

Although the prediction of protein localization from the amino acid sequence is beyond the scope of this overview, some published methods are worth mentioning. A mainly empirical method originating from experimental and computational observations, called PSORT, has been developed since the beginning of the 1990s [99] and is available at http://psort.ims.u-tokyo.ac.jp. More mechanistic approaches, i.e. simulating the mechanism of cellular sorting, are MultiLoc [56] (available at (http://www-bs.informatik.uni-tuebingen.de/Services/MultiLoc), LOCtree [98] (available at http://cubic.bioc.columbia.edu/services/loctree), and PLOC [105] (available at http://www.genome.jp/SIT/plocdir).

8 Exchange of Model Data

In the last two decades, with the rapidly increasing number and complexity of publicly available models, efforts have been made to introduce standards for the flexible exchange and re-use of biological model descriptions between different simulation packages. These efforts try to minimize the need for re-formulation of model implementations when switching between different simulation and/or analysis tools. They involve both the development of multi-purpose unified model description languages (described below), as well as the need for the definition of standards for the model architecture and implementation, such as the necessity for a single reference description, the traceability of the choice of parameters, initial conditions etc. [60]. A widely accepted standard is Minimum information requested in the annotation of biochemical models (MIRIAM) [102]. Additionally, published model implementations in public databases such as CellML are often subject to increased requests for annotations and modification history documentation. Also, equations and identifiers are unified or re-annotated if possible [81] according to ontological frameworks [140].

The following sections will focus on the most important standard model exchange formats. For further overviews of this topic, the reader is referred to dedicated reviews [125]. A variety of modeling software packages will be summarized in Sect. 9.

8.1 SBML

The systems biology markup language (SBML) was introduced in 2001 and extended to SBML level 2 in 2002 [58], the most recent version 4 having been published in 2008 [59]. It is based on XML (extensible markup language), a widespread general hierarchical document description format. It provides a generalized description of biochemical reaction systems, e.g. cell signaling pathways, metabolic pathways, biochemical reactions, gene regulation etc.

Since SBML is strictly standardized and commonly used, there is a variety of helpful applications available in the web, e.g. an automated SBML validator that checks for syntax and consistency (http://sbml.org/Facilities/Validator), a converter into human-readable documentation formats (PDF, LaTeX etc., http://www.ra.cs.uni-tuebingen.de/software/SBML2LaTeX/index.htm) [28], automated creation of kinetic equations (http://www.ra.cs.uni-tuebingen.de/software/SBMLsqueezer) [27] and others.

An increasing number (over 450 by August 2010) of peer-reviewed, quantitative models of biochemical and cellular systems is collected in the BioModels database (http://www.ebi.ac.uk/biomodels-main) [101]. To be accepted in BioModels, the model must comply with MIRIAM [102]; additionally, the model components are annotated according to standard databases (UniProt, KEGG, Reactome etc).

8.2 CellML

Another popular model description language is CellML [80], which has been developed starting in 1998. The first specification was published in 2001 [49]. Like SBML, it is XML-based; however, it is more general in that it is not restricted to biological systems.

The CellML repository (http://www.cellml.org/models) currently holds over 480 models (as of August 2010), mostly derived from peer-reviewed literature, covering a wide range of biological processes such as signal transduction pathways, metabolic pathways, electrophysiology, immunology, cell cycle, muscle contraction, mechanical models and constitutive laws etc. [81].

As with SBML, a variety of helper tools is available for CellML, e.g. converters to SBML: CellML2SBML [117] etc. A review of CellML-associated software can be found in [37].

9 Cell Modeling Program Packages

A variety of simulation software packages is available for the facilitation of whole-cell simulations or for solving specific biological questions without the need for extensive programming and mathematical and numerical optimization efforts. This overview concentrates on six different packages that are applicable to simulation of mesoscale intracellular transport processes: VirtualCell, MCell, Bio-SPICE, StochSim, COPASI and BioDrive. Not included in this overview are simulators for multicellular systems [55], pure genomic network-based simulation environments such as ECell [131]/ ECell2 [127] and CellDesigner [35], or tools that are dedicated for neural simulation, e.g. NEURON [54] or GENESIS [8]. Also, simulation suites for molecular dynamics simulations, such as CHARMM [9], AMBER [14] and GROMACS [53], which are to some extent suitable for simulation of transport processes at the nanoscale level [5, 68], are beyond the scope of this overview.

9.1 Virtual Cell

The Virtual Cell, or VCell, project, which was first published in [116], is developed at the National Institutes of Health (NIH). It is available from http://vcell.org. It provides a framework for modeling biochemical, electrophysiological and transport phenomena and allows the user to build complex models with a Java-based interface. Arbitrary (but temporally constant) compartment topologies and geometries are supported. The computation is performed on-line on computing clusters dedicated to VCell at the National Resource for Cell Analysis and Modeling (NRCAM).

VCell biological models are composed of three components—a physiological model containing the mechanistic hypothesis (i.e. reactions, fluxes, electrical currents on membranes), an application with experimental conditions, geometry and modeling approximations, and one or more simulations (Fig. 5). The underlying mathematics and physics as well as the working principles are described in [94, 121]. For improved visualization of simulation results on physiological data, a 3D interactive visualization tool is integrated.

Fig. 5
figure 5

General VCell workflow: from general physiology, over specific structures and applications, to simulation. a Physiology: topology and reactions/fluxes. b Application: structure mapping and boundary conditions. c Simulation

Full size image

9.2 MCell

MCell, available at http://www.mcell.cnl.salk.edu, is a general Monte Carlo simulator of cellular microphysiology. It enables simulation of (intra-)cellular signalling in and around cells. The framework utilizes Monte Carlo random walk and chemical reaction algorithms based on pseudo-random-number generation and is thus capable of tracking the stochastic behavior of discrete molecules in space. A prominent feature of the software is its ability to stop the simulation at arbitrary time positions, change parameters and morphologies on-line and continue with the simulation, which in principle allows for modeling of dynamic cell shape variations. The software utilizes a special model description language (MDL) and currently does not provide an interactive interface.

Although originally having been developed for miniature endplate current generation in the vertebrate neuromuscular junction [7], this software suite is meanwhile generalized enough to be used for different biological questions, such as \(\hbox{Ca}^{2+}\) dynamics in dendritic spines [33] and modulation of impulse propagation in nodes of Ranvier depending on \(\hbox{Na}^+\)channel distribution [124]. The simulations can be performed on distributed computing environments—the “grid” [32]—to increase computation speed [13].

9.3 Bio-SPICE

Bio-SPICE is an open-source project (http://biospice.sourceforge.net), originally developed at the Defense Advanced Research Projects Agency (DARPA) in 2002. It is intended for modeling and simulation of spatiotemporal processes in living cells. The toolbox combines different software from various vendors by integrating them into the “Dashboard”, the Bio-SPICE core application that provides a consistent workflow for modeling, analysis and simulation. It is closely related to the Systems Biology Workbench (SBW), a software framework that allows heterogeneous application components written in multiple programming languages on different platforms [115].

9.4 StochSim

StochSim, a discrete stochastic simulator for chemical and biochemical reactions [31], available at http://www.pdn.cam.ac.uk/groups/comp-cell/StochSim.html, with a graphical user interface at http://www.ebi.ac.uk/lenov/stochsim.html, has been developed as part of a study of bacterial chemotaxis [95]. This algorithm considers chemical reaction partners as individual interacting objects. The StochSim algorithm is in most cases more efficient than the older Gillespie [38, 39] algorithm, especially when multi-state molecules are considered [92].

9.5 GEPASI and COPASI

For simulation and analysis of biochemical reaction networks with support of recent model exchange standards such as SBML and CellML (Sect. 8), the COmplex PAthway SImulator (COPASI) software suite has been developed [57, 91] as a successor to the GEneral PAthway SImulator (GEPASI) [90] that was first published in the early 1990s. It is available for multiple platforms underhttp://www.copasi.org. The software considers sensitivity and metabolic control analysis methods as well as various optimization algorithms such as evolutionary programming, Nelder–Mead, particle swarm, simulated annealing and others [91].

9.6 BioDrive

The BioDrive biochemical reaction and gene expression simulation software suite has been developed by the group of Kyoda et al. since the late 1990s. The authors emphasize its easy usability, or “biologist-friendliness”, as well as its ability to cope with multi-cellular organisms and extra-cellular processes. It is based on ordinary differential equations and incorporates diffusion and spatiotemporal patterning [75].

10 Conclusions

In the past two decades, significant efforts have been made to deliver more in-depth insights into the dynamics of mesoscale intracellular transport processes in eukaryotic cells and their impact on cell function. A large variety of general-use biochemical simulation software is available, though only a subset (VCell, MCell, Bio-SPICE) considers spatially resolved models. Active transport via, e.g., the actin filament or microtubule network is not yet explicitly incorporated. The influence of such processes has increasingly moved into the scope of research, hence numerous more or less ad-hoc Monte Carlo methods have been implemented and published. Unfortunately, the interchangeability of such models is limited compared to standardized model description languages. The unification of different models covering specialized aspects to more holistic approaches with a perspective towards whole-cell simulations is therefore still a very cumbersome process; however, the foundations for the necessary infrastructure are laid (e.g. MIRIAM) and are being further extended.