Agent-Based Modeling of Intracellular Transport
MMirko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport
European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255
Agent-Based Modeling of Intracellular Transport
Mirko Birbaumer, Frank Schweitzer
Chair of Systems Design, ETH Zurich, Kreuzplatz 5, 8032 Zurich, Switzerland
Dedicated to Werner Ebeling on the occasion of his 75th birthday
Abstract
We develop an agent-based model of the motion and pattern formation of vesicles. Theseintracellular particles can be found in four different modes of (undirected and directed)motion and can fuse with other vesicles. While the size of vesicles follows a log-normaldistribution that changes over time due to fusion processes, their spatial distribution givesrise to distinct patterns. Their occurrence depends on the concentration of proteins which aresynthesized based on the transcriptional activities of some genes. Hence, differences in thesespatio-temporal vesicle patterns allow indirect conclusions about the (unknown) impact ofthese genes.By means of agent-based computer simulations we are able to reproduce such patterns onreal temporal and spatial scales. Our modeling approach is based on Brownian agents with aninternal degree of freedom, θ , that represents the different modes of motion. Conditions insidethe cell are modeled by an effective potential that differs for agents dependent on their value θ . Agent’s motion in this effective potential is modeled by an overdampted Langevin equation,changes of θ are modeled as stochastic transitions with values obtained from experiments,and fusion events are modeled as space-dependent stochastic transitions. Our results for thespatio-temporal vesicle patterns can be used for a statistical comparison with experiments.We also derive hypotheses of how the silencing of some genes may affect the intracellulartransport, and point to generalizations of the model. Agent-based modeling has proven to be a versatile tool to simulate processes of structure for-mation bottom up. By assuming features and interaction rules of agents on the “microscopic"level, one is able to observe the emergent systems properties on the macroscopic level. This isof particular importance in those areas where the systems dynamics can hardly be captured topdown, i.e. in living systems, including biological, social or economic systems.But the advantage of agent-based models in freely defining agent properties and interactions soonturns out to be a pitfall, because this way arbitrary patterns can be generated and it is difficultto choose the right values in a high dimensional parameter space. To minimize these problems,there are basically two ways: (i) to closely link the agent’s properties to experimentally observeddata, and (ii) to apply methods that allow to aggregate the agent dynamics, to formally derive1/21 a r X i v : . [ q - b i o . S C ] S e p irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 the systems dynamics. The latter provides a firm relation between agent’s features and systemsfeature’s and may reveal also the role of certain (control) parameters.The concept of Brownian agents [16] was developed to facilitate the second way. The dynamics ofagents is described in a stochastic manner, similar to the Langevin approach of Brownian motion.This allows to obtain on the macroscopic level a closed-form partial differential equation for thedensity, that for the case of Brownian motion simply describes a diffusion process. The dynamicsin most real systems, however, is much more complicated. Agents are not simple random walkers,they respond to information in their environment, follow chemical gradients, and can at the sametime also contribute to generating information, chemical gradients etc. Further, agents do notbehave the same all the time. Instead, they may have different modes of activity each of whichcorresponds to a particular behavior. To cope with these features, Brownian agents are describedby internal degrees of freedom and their environment is modeled as an adaptive landscape, oreffective potential, which can be modified by the agents while responding to the informationprovided. Further, transitions between the agent’s internal degrees are possible, dependent oninternal or external conditions.On the formal level, the macroscopic dynamics is then no longer described by a diffusion equation,but by a quite advanced reaction-diffusion equation with a variable drift term, which is coupledto another differential equation describing the dynamics of the adaptive landscape dependent onthe agent’s activity. This allows to tackle the dynamics of systems comprised of many interactingagents on two levels: (i) the agent level, where computationally efficient computer simulations canbe performed, (ii) the system’s level, where coupled differential equations may be obtained andinvestigated analytically. The application discussed in this paper is unfortunately complex enoughto not provide closed form equations for an analytical treatment. Nevertheless, the concept ofBrownian agents allows us to formally specify the agent dynamics in terms of stochastic equationsof motion in an adaptive landscape.The aim of this paper is twofold: (i) to develop an agent-based model of intracellular transport andpattern formation, which is general enough to be applied to various such phenomena involvingfree and directed motion and fusion processes (Section 3), and (ii) to specify this agent-basedmodel for the case of vesicle movement and fusion, in close relation to experimental findings(Section 4). Importantly, the internal degrees of freedom of agents and transitions between theseare obtained from experiments. This allows us to observe pattern formation on real time andspatial scales (Section 5), the outcome of which can, at least in a statistical manner, be comparedwith real experiments. Hence, applying the concept of Brownian agents to a real problem, i.e.the intracellular transport and pattern formation of vesicles, demonstrates both the versatilityof the concept and its suitability to generate hypotheses about real intracellular processes.2/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport
European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255
In this paper, we are interested in the intracellular transport and pattern formation of vesicles .These are quite small intracellular particles (diameter approx. . µ m) (see [15], [2] and [6]).They are formed at the cell membrane, to contain some extracellular material engulfed by thecell membrane. This import of material, called endocytic cargo, may include macromolecules,but also viruses, which are all encapsulated in vesicles – a process called endocytosis (see [6], [13]and [14]). In this paper, we do not consider endocytosis explicitly, but assume that vesicles areformed at the membrane and then released into the interior of the cell at a constant rate (latercalled internalization rate). It is known from experiments that the size distribution of newlyformed vesicles follows approximately a log-normal distribution (see [1]).Figure 1: Two-dimensional representation of a cell with cytoskeletal structure. Adapted from[4].Released vesicles can diffuse inside the cell, but they can also be reabsorbed by the membraneat a constant rate, later called recycling rate ([3]). Vesicles need to be transported from the cellmembrane to the endosomal system located in different areas inside the cell – this transportprocess is called membrane trafficking (see [6]). To facilitate this process, in addition to the freediffusion , vesicles can also perform a directed motion along the cytoskeleton , which is an intracel-lular structure made of two different kinds of polymer filaments: actin filaments and microtubulefilaments (see Figure 1). Whereas actin filaments are randomly distributed, microtubules (MT)are directed toward the microtubule organization center (MTOC), which is located close to thecell nucleus (see Figure 2). In order to perform a directed motion along actin or MT filaments,vesicles have to bind to motor proteins (kinesines, dyneins and myosins) which basically deter-mine the type of motion. The interaction between motor proteins and their related filaments3/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 depends on several additional factors, most notably the presence of specific proteins such as RabGTPases, scaffolding proteins and receptor proteins (see [9]). To which extent these proteins arepresent depends on the other hand on the genes of the cell, which have to be transcriptionallyactive in order to synthesize these proteins – a process called gene expression .It is exactly this dependency, which motivates our interest in the motion of vesicles. If a givengene, for example cdk8 , is transcriptionally active, it allows the synthesis of the protein CDK8,which may affect intracellular processes in various, mostly unknown ways. In this paper, weprecisely ask how such a gene – through the synthesis of the specific protein – affects the motionand pattern formation of vesicles. The latter process results from the fact that vesicles can formlarger ones by fusing with other vesicles. The fusion process relies on energy provided by the celland can only take place on the MT if vesicles are close enough and below a critical size. Thetwo simultaneous dynamic processes, namely (free or directed) motion and fusion result in adistinct spatio-temporal distribution of vesicles of different sizes - which we want to predict withour model.The patterns produced by means of our agent-based simulation can then be compared to experi-ments which show such vesicle patterns depending on the transcriptional activity of specific genes(see [1]). These genes can, for example, be knocked out in RNAi or drug screens, which in turnperturb the synthesis of proteins (see [11]). Of course, such patterns can be only compared in astatistical sense, a problem discussed in the Conclusions. However, if we are able to reproduceempirical patterns with our model, we argue that the underlying dynamic processes, motion andfusion, are covered sufficiently with our modeling assumptions. This does not only hold for theassumed interaction between vesicles and MT or other vesicles, it shall also hold for particularparameter dependencies, most notably the concentration of specific proteins. Precisely, we wantto end up with a testable prediction of how these concentrations affect vesicle motion and fusion,which shall be confirmed by subsequent experiments (see [1]). Some of the transition rates laterused in our model specifically depend on the experimental setup, e.g. the endocytic cargo. Herewe consider the case of Transferrin, an iron-binding protein contained in the vesicles to which flu-orescently labelled proteins can be attached, i.e., vesicles can be made visible in the experiment.The aforementioned internalization rate and recycling rate are thus taken for Transferrin.Eventually, we note that our modeling approach (as every modeling approach) is based on anumber of simplifications: (i) we neglect the motion of vesicles along actin filaments, because itwas shown [4] that such processes do not affect the pattern formation (recall that fusion onlytakes place on MT), (ii) we assume that the cytoskeleton is described by the spatial structureof MT only (i.e. tha actin filaments are neglected) and that MT are abundant (i.e. there arealways MT to move on) and do not change in time. This allows to neglect the growth andshrinkage of MT when modeling the motion and fusion of vesicles. (iii) We neglect fission, i.e.the fragmentation of larger vesicles into vesicles of smaller sizes.4/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport
European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255
Brownian agents
Our modeling approach is based on the concept of Brownian agents [16]which found many interesting applications in biology [7, 8, 10, 12] but also in modeling socialsystems such as online communities [17]. It allows to formalize the agent dynamics using methodsestablished in statistical physics. A Brownian agent is described by a set of state variables u ( k ) i ,where the index i = 1 , ..., N refers to the individual agent i , while k indicates the different vari-ables. These could be either external variables that can be observed from the outside, or internaldegrees of freedom that can only be indirectly concluded from observable actions. Noteworthy,the different (external or internal) state variables can change in the course of time, either due toinfluences of the environment, or due to an internal dynamics.In the following, each agent represents an intracellular vesicle which, in accordance with theprevious description, is able to change its state by spatial mobility, changes of activity, andgrowth processes, as formalized subsequently. Spatial mobility
For the agent’s spatial position r i ( t ) , We assume that changes in the courseof time can be described by an overdamped Langevin equation of a Brownian particle moving inan effective potential [16]: d r i dt = v i [ θ i ( t )] = α [ θ i ( t )] γ ∂h e ( r , t, θ ) ∂ r + √ D ξ i ( t ) (1)The overdamped limit implies that the absolute value of the agent’s velocity v i is approximatelyconstant, but the direction may change due to stochastic influences. Further, v i implicitly dependson the agent’s mode of activity, θ i ( t ) .Eqn. (1) assumes that the agent’s motion is influenced by two different forces, a deterministic onewhich results from the gradient of the effective potential, and a stochastic one, which is assumedto be Gaussian white noise, (cid:104) ξ i ( t ) (cid:105) = 0 , (cid:104) ξ i ( t ) ξ j ( t (cid:48) ) (cid:105) = δ ij δ ( t − t (cid:48) ) . The strength of the stochasticforce D = k B T /γ determines, in the spatial case, the diffusion coefficient, with γ being thefriction coefficient.The deterministic part contains two important ingredients: The effective potential h e ( r , t, θ ) de-scribes the conditions inside the cell. The response function α [ θ i ( t )] depends on the internal stateof the agent, θ i ( t ) and determines what component of the effective potential actually influencesthe agent. Both are specified later after we made clear the notion of the internal state θ . Changes of activity
We assume that the agent’s mode of activity θ i ( t ) can be changed, ω ( θ (cid:48) | θ ) being the transition rate from state θ into any other state θ (cid:48) . In accordance with the literature5/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 θ = 0 Free diffusion in the cytosol θ = 1 Kinesin-driven transport towards MT plus-ends θ = 2 Dynein-driven transport towards MT minus-ends θ = 3 Kinesin/Dynein-driven transport on MT with tendency to fusion with other vesicles θ = 4 Bound to the cell membraneTable 1: Different modes of activities, MT stands for microtubuli[4], we distinguish five different modes of activity of a vesicle as shown in Table 1, each of whichis expressed by a different value of the internal degree of freedom θ .While in principle transitions between all states could be assumed, only a subset of them isbiologically relevant. Table 3 in Section 4.1 will list those together with their respective value,i.e. the expected number of transitions per time unit. Growth and decay
We assume that fusion, i.e. the coalescence of two vesicles with sizes s i and s j , can be described by a transition rate ω ( s i + j | s i , s j ) that depends on the internal states θ j , θ i of the agents, i.e. their ability to fuse, and their effective distance | r i − r j | . As describedabove, fusion is only possible for agents with the internal state θ = 4 which is expressed by theKronecker delta, δ θ i , . Further, due to the volume exclusion, agents can effectively approach eachother only up to a distance d (which represents an average spatial extension of vesicles). Theability to fuse also depends on the vesicle size because of the energy required for this process.Because the available fusion energy is limited, it was observed experimentally that vesicles ofa size larger than s max do not fuse. This is considered in the transition rates by an additionalexponential cut-off term which becomes effectively zero if one of the fusing vesicles reaches themaximum size. This leads us to the transition rate for fusion: ω ( s i + s j | s i , s j ) = ω s δ θ i , δ θ j , d + | r i − r j | · e ε (2 s max − s i − s j ) (cid:2) e ε ( s max − s i ) (cid:3) · (cid:2) e ε ( s max − s j ) (cid:3) (2)Here ω s denotes the fusion affinity which depends on the protein concentration, and ε = 0 . ischosen as a small number, to increase the cutoff effect.Little is known regarding the fission process, i.e. the fragmentation of a vesicle of size s l into twovesicles of sizes s i , s j (with s l = s i + s j ) . Therefore, we assume that the respective transition rate ω ( s i , s j | s i + s j ) is a constant w equal for all possible fragmentation processes, which describesspontaneous fragmentation. If w is small compared to other transition rates, fission can beneglected in first approximation.We note that, because of the fusion process, the total number of vesicles, is no longer constant.While a conservation of the total mass M of all vesicles can still be assumed, both the number of6/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 vesicles and their size distribution P ( N , N , ..., N s ..., t ) changes over time. Therefore, we have M = N (cid:88) l =1 s N s ( t ) = const. ; N (cid:88) l =1 N s ( t ) (cid:54) = const. (3) Effective potential
After the above distinction between the different values for the internaldegree of freedom θ , we can now specify the effective potential which depends on these states. h e ( r , t, θ ) denotes a scalar potential field that results from the influence of different field compo-nents h θ ( r , t ) . Each of these components refers to specific conditions inside the cell. Comparedto the time scale involved in the motion of the vesicles, some of these conditions can be assumedas constant in time, but varying across space.With reference to Table 1, h ( r ) describes the influence of the cell membrane on the motion ofagents in state θ = 4 as they can bind to the membrane. h ( r ) and h ( r ) determine the agent’smotion along the microtubule filaments. h ( r ) on the other hand represents the cell topology,i.e. it generates a repelling force close to the cell membrane and to the nucleus, but is is simplya constant inside the cell, because free diffusion inside the cell should not be affected.The only time-dependent component of the effective field is h (∆ r , t ) , which affects the fusionprocesses between vesicles (fission neglected). In fact, this is a short range attraction potentialwhich increases with decreasing distance ∆ r = | r i − r j | . Since agents move, h changes in timedepending on their actual positions r i ( t ) and internal states θ i ( t ) .In order to describe how the effective potential results from the different field components, wehave to consider the response function α [ θ i ( t )] that determines which of the field components areactually “seen” by the agents conditional on their internal states. In accordance with the abovedistinction, we specify: α [ θ i ( t )] h e ( r , θ, t ) = α h ( r ) + δ θ i , α h ( r ) + δ θ i , α h ( r )+ δ θ i , α h ( r ) + δ θ i , α h (∆ r , t ) (4)The different α k are dimensional constants. With eqn. (4) the motion of every agent is specifiedaccording to eqn. (1). We point out that agents with θ i = 0 behave like simple Brownian particleswhich are not affected by any conditions inside the cell, except for the boundary conditions. Master equation
The state of each individual agent is now described by a triple of threedifferent state variables { r i ( t ) , θ i ( t ) , s i ( t ) } that can change in time according to the processesspecified above. The multi-agent system is thus described by the grand-canonical N -particledistribution function P N = P N ( r , θ, s, t ) = P N ( r , θ , s , ..., r N , θ N , s N , t ) , (5)7/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 which describes the probability density of finding N Brownian agents with the distribution ofinternal parameters θ , positions r and sizes s at time t . Note that N is not constant but canchange over time due to fusion (and fission) events.The complete dynamics for the ensemble of agents can be formulated in terms of a multivariatemaster equation : ∂∂t P N ( r , θ, s, t ) = − N (cid:88) i =1 (cid:26) ∇ i (cid:20) α [ θ i ( t )] γ ∇ i h e ( r , t, θ ) P N ( r , θ, s, t ) (cid:21) − D ∆ i P N ( r , θ, s, t ) (cid:27) (6a) + N (cid:88) i =1 (cid:88) θ (cid:48) i (cid:54) = θ i (cid:2) ω ( θ i | θ (cid:48) i ) P N ( r , θ (cid:48) i , θ ∗ , s, t ) − ω ( θ (cid:48) i | θ i ) P N ( r , θ, s, t ) (cid:3) (6b) + N (cid:88) i =1 (cid:88) i According to the distinction between the activity modes given in Table 1, the movement ofthe agents can be of two types: free motion by diffusion processes, described by the diffusioncoefficient D , and bound motion along the microtubuli. Both types are described by eqn. (1). Inorder to calculate the different v i ( θ ) , we would need to specify explicitly the related componentsof the effective potential, h θ ( r , t ) that result in the directed motion along the microtubuli. Tosimplify the procedure and match it with experimental findings, we may instead consider thevalue of v i ( θ ) as given by experiments. Then, the role of the field component is reduced to simplykeep the agents on the microtubuli as long as they are in the respective internal states θ . I.e.both ∇ i h ( r ) and ∇ i h ( r ) just determine the direction of motion along the microtubulus towhich agent i is attached, whereas h ( r ) specifies the boundary condition given by the locationof the cell nucleus and the cell membrane (see also Figure 2). Table 2 provides the values for thevelocities which are assumed to be constant and the same for all agents in the respective internalstate. State Parameter value Description θ = 0 D = 10 − m / s Diffusion coefficient in cytoplasma θ = 1 v = 0 . µ m / s Velocity of vesicle on MT towards plus-ends θ = 2 v = 0 . µ m / s Velocity of vesicle on MT towards minus-endsTable 2: Parameters describing the free ( θ = 0 ) or bound ( θ = 1 , ) motion of agents. MTmeans ‘microtubulus’. Note that the diffusion coefficient is related to the friction coefficient via γ = k B T /D , which results in γ = 4 · − kg / s . Values according to [5].As a next step, we need to specify the transition rates ω ( θ (cid:48) | θ ) between different modes of activ-ity. Again, instead of providing explicit expressions, we may simply take values obtained fromexperiments. These values are available to us only for the unperturbed state of the cell, i.e. for abaseline or reference case in which conditions inside the cell are not changed on purpose. In thefollowing, baseline values are indicated by the superscript b . Table 3 lists all biologically relevanttransition rates for the unperturbed scenario, as obtained either from the literature or from ownexperiments. The given set of transition rates leaves us with a large degree of freedom which however turns outto be a pitfall: if we wanted to compare the simulations with patterns observed in the experiment,9/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 Parameter Value Description ω b (1 | =0.05 s − Transition from diffusive to MT (plus-end) bound state ω b (0 | =0.33 s − Transition from MT-bound (plus-end) to diffusive state ω b (2 | =0.05 s − Transition from diffusive to MT (minus-end) bound state ω b (0 | =0.33 s − Transition from MT (minus-end) bound to diffusive state ω b (3 | =0.01 s − (cid:63) Transition from MT-bound to MT-bound state with fusion ω b (1 | =0.02 s − (cid:63) Transition from MT-bound state with fusion to MT-binding ω b (3 | =0.01 s − (cid:63) Transition from MT-bound to MT-bound state with fusion ω b (2 | =0.02 s − (cid:63) Transition from MT-bound state with fusion to MT-binding ω b (0 | =0.0025-0.0033 s − Internalization rate of vesicles ω b (4 | =0.00083-0.0012 s − Recycling rate of vesiclesTable 3: Biologically relevant transition rates ω b ( θ (cid:48) | θ ) for the baseline case (unperturbed con-centrations of proteins). Values according to [5], values with (cid:63) are from own experiments, see[1].we would need to raster the full parameter space to find appropriate combinations of transitionrates which lead to a realistic outcome. There are basically two ways to reduce this parameterspace: (i) we use known values of the transition rates as e.g. reported in the literature, (ii) weintroduce reduced transition rates, assuming that not all processes are really independent. Wefollow a combination of the two, defining the reduced transition rates as given in Table 4 Ω = ω (1 | ω (0 | Affinity for microbule plus-ends Ω = ω (2 | ω (0 | Affinity for microbule minus-ends Ω = ω (3 | ω (1 | Fusion tendency on microbule plus-ends Ω = ω (3 | ω (2 | Fusion tendency on microbule minus-ends Ω = ω (0 | ω (4 | Internalization versus recycling rateTable 4: Reduced transition rates Ω i , to combine two transition rates ω which refer to thesame intracellular (transport) mechanism but describe inverse processes.In the following we assume that Ω = Ω ≡ Ω F which denotes the fusion tendency that issupposed to be independent of the kind of motor protein involved (this was kinesin for Ω and10/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 dynein for Ω ). With this, we finally have reduced the transition rates to 4 parameters given by Ω , Ω , Ω F and Ω .These reduced parameters of course do not fully determine the individual transition rates whichare needed to recover the correct dynamics. Therefore, for one of the transition rates we use thebaseline value from the literature, as given by Table 3. So, we define e.g. Ω = ω (1 | /ω b (0 | or Ω = ω b (1 | /ω (0 | .Eventually, we want to emphasize the distinction between the reduced transtion rates Ω ui forthe unperturbed case (control cell) and Ω ki for the perturbed case, where the protein k wasmanipulated. Similar to the discussion in Section 4.1, they refer to the same transition rates ω ,but with different concentrations c k . In order to complete the setup for the stochastic computer simulations, we still need to specifythe boundary conditions which refer to the cell geometry. Figure 2 (a,b) shows the two differentcell geometries chosen, a rather circular cell and an elongated one. The outer boundary of thecell membrane and the inner boundary of the nucleus are both assumed to be impermeable walls,described by a hard sphere potential h ( r ) .The interior of the cell contains the cytoskeleton, i.e. both actin filaments and microtubuli (MT)which provide boundary conditions for the motion of vesicles (see also Figure 1). As alreadydiscussed, we do not consider motion along actin filaments because the related processes do notcontribute to vesicle pattern formation. MT, on the other hand, are abundant and always pointto the microtubule organization center (MTOC) which is assumed to be in the perinuclear areaon the side with the largest part of the cytosol (see Figure 2 a,b).Regarding the initial conditions, we need to specify the number , the position , the internal state and the size of the agents at t = 0 . For our model, we assume that initially all vesicles are boundto the cell membrane, i.e. at t = 0 agents start with θ i (0) = 4 at a position r i (0) randomlychosen from the cell boundary. In order to determine the initial number of agents, we startfrom the experimental observation that an average cell (of the type considered) contains about200 internalized vesicles at steady state. This number excludes vesicles still bound to the cellmembrane. We further know from experiments that internalization rates, i.e., transition ratesfrom the initial bound into the free moving state (and vice versa), i.e. ω (0 | 4) = 3 . · − s − and ω (4 | 0) = 8 . · − s − (see Table 3). If N (0) = N (0) denotes the (unknown) number of agentsat t = 0 (all assumed to be in the bound state) and N (0) − N st = 200 the (known) number ofagents no longer in the bound state at steady state ( st : steady state), then we can postulate for11/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 MTOC MTOC Figure 2: (left) (a) Circular geometry of a non-treated (control) cell; (right) (b) Elongated ge-ometry of a perturbed cell (treated with siRNA, which leads to the depletion of CDK8). Thescale bar corresponds to 5 µm . The microtubule filaments are schematically sketched, pointingtoward the MTOC (microtubule organization center). Vesicle in an internal state θ ∈ { , , } move along microtubules in both directions.the dynamics for N ( t ) the following rate equation: dN ( t ) dt = − ω (0 | N ( t ) + ω (4 | 0) [ N (0) − N ( t )] (8)After a time t ∼ / [ ω (0 | 4) + ω (4 | , this dynamics reaches a steady-state solution N st = ω (4 | ω (0 | 4) + ω (4 | N (0) (9)from which we can calculate N (0) assuming that N (0) − N st = 200 . This approximation neglectsthe fact that the total number of vesicles have decreased at time t because of the fusion betweenvesicles. Thus, we may slightly increase the initial number of agents and have eventually chosen N (0) = 350 .It remains to determine the initial distribution of vesicle sizes. We want to start from a mostrealistic one because (a) later we want to compare the time scale of structure formation with theexperimental observation, and (b) because our modeling setup has neglected the fragmentationrate of vesicles (which would be needed if an arbitrary initial distribution needs to relax into12/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 a realistic one). Again, we rely on experimental observations [1] that have found a log-normaldistribution of vesicle sizes: P ( s ; µ, σ ) = 1 sσ √ π exp (cid:26) − (ln s − µ ) σ (cid:27) , (10)where µ and σ are the mean and standard deviation of the variable’s natural logarithm. We now have determined all ingredients for stochastic computer simulations which include thefollowing dynamic processes (specified on the agent level): Initialization At t = t , 350 agents with θ i (0) = 4 are randomly placed at the cell boundary.Their initial size s i (0) is drawn from the log-normal distribution, (10). Movement Agents can change from the bound state into the free moving state at a rate ω (0 | and from the free moving state into directed motion at rates ω (1 | , ω (2 | . They all moveaccording to the equation of motion, (1).In our modeling approach, only agents with the internal states θ ∈ { , , } move alongthe MT. We assume that, whenever an agent switches from θ = 0 into either θ = 1 or θ = 2 , i.e. from a free motion into a bound motion, a MT is “always at hand”. If the agentswitches from θ ∈ { , } into θ = 3 where it is ready to fuse, it continues to move into thesame direction as before, until it collides with another agent. Fusion Precisely, fusion occurs only on MT. Agents need to be in state θ = 3 and in a sufficientlyclose distance. After fusion, the smaller agent “disappears”, whereas the larger one hasincreased its size, but keeps the internal state θ = 3 until one of the transitions ω (2 | or ω (1 | happen. Inversion Agents bound to the MT can become freely moving at the rates ω (0 | , ω (0 | ,whereas free moving agents can be bound again to the cell membrane at the rate ω (4 | . Concurrency We assume that agents can move and transit into different internal states at thesame time. This allows us to decouple the motion of agents from the various reactions(change of internal states and fusion).The state of the multi-agent system is at any time completely described by the N -particledistribution function P N ( r , θ, s, t ) . However, because of the dynamical processes, the systemstate always changes and its average “life time” T m is just the inverse of the sum over all possible13/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 transition rates that can change the given state (including changes of position, internal states,and sizes).Because we have to deal with movement and reactions at the same time, we have chosen asufficiently small fixed time interval ∆ t = 0 . s to solve the equations of motion of each individualagent. To answer the question if in the respective time interval also a change of the internal statesor a fusion process occurs, we proceed as follows:Each of the possible reactions has a different probability to occur, which is determined by therespective transition rate and the available time, ∆ t . In order to pick one of the possible reactions,we draw a uniformly distributed random number U ∈ RND [0 , and choose the process z thatsatisfies the condition z (cid:88) j =0 ω z ( · )∆ t < U < z +1 (cid:88) j =0 ω z ( · )∆ t (11)Because ∆ t was chosen such that (cid:80) Nn =1 ω n ( · )∆ t ≤ , none of the possible processes is excludedfrom being picked. On the other hand, it may occur that none of the processes is being chosen ifthe sum is much smaller than 1 and U close to 1. Then no reaction occurs during the respectivetime interval, but movements take place. Figure 3 presents computer simulations of the vesicle patters for both the unperturbed (control)cell and the perturbed cell (cf also Figure 2). We emphasise that these simulations refer to realspatial and time scales, so they should be comparable, at least in a statistical sense, to patternsobserved from experiments. These experiments are reported in [1] and have motivated the choiceof the reduced transition rates, Ω i , which are treated in this paper as free parameters.Comparing the pattern formation in the perturbed cell with the one in the control cell, we note anumber of differences: we observe a localization of large vesicles on the one hand at the tips of theelongated branched-out perturbed cell, on the other hand large vesicles are as well localized inthe perinuclear area of this cell. In the unperturbed cell, the majority of large vesicles are locatedaround the nucleus and vesicles are spread over the entire cell surface getting more sparse towardsthe cell periphery.Comparing the vesicle size distribution of both the perturbed and the unperturbed cell, we findthat they preserve the form of the log-normal distribution, but the mean value µ , in the courseof time, shifts to significant larger values in the perturbed cell (see also Table 5). The perturbedcell displays a geometry that may have facilitated fusion of vesicles in its branches within which14/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 they accumulate. The transport of vesicles from these branches to the nucleus of the perturbedcell seems to be suppressed. t µ u σ u µ σ µ ( t ) and standard deviation σ ( t ) of the log-normal distribution at differenttimes t (min) for the unperturbed (superscript u ) and the perturbed cell. As pointed out above, we found very realistic vesicle patterns both for the perturbed and thecontrol cell for the given set of reduced transition rates Ω i (see Figure 3). These transition ratesare treated as free parameters in our model – but provided they are correct (which can onlybe confirmed by comparing the patterns with experiments, statistically) they allow an indirectestimation of the concentration dependence of the transition rates ω ( θ (cid:48) | θ ) .As already stated, the transition rates are given to us only for the unperturbed state of the cell.Table 3 presents the values for the baseline case. One should note, however, that the baseline valuenot necessarily describes the experimental situation because it was obtained under conditionswhich are hardly reproducible completely. If we for example change the concentration c k of someprotein k inside the cell which is involved in processes of fusion or directed motion, this willcertainly change the value of the respective transition rates, i.e. ω = ω ( c k , c l , ... ) . Hence, it is notonly sufficient to know the values of the baseline case, we should also know how these values ofthe transition rates change with the concentrations c k . If we denote the transition rates in theunperturbed case by ω u (omitting the θ dependence at the moment) and the respective proteinconcentrations by c u we may assume in first-order approximation the following expansion: ω ( c k , c l (cid:54) = k = c ul (cid:54) = k ) = ω u ( c uk , c ul (cid:54) = k ) + ∂ω∂c k (cid:12)(cid:12)(cid:12)(cid:12) c uk ( c k − c uk ) . (12)To further specify the functional dependency ω ( c k ) we make the following ansatz: ω ( c k , c l (cid:54) = k = c ul (cid:54) = k ) = ω u · exp (cid:18) κ k c k − c uk c uk (cid:19) . (13)which satisfies ω ( c k , c l (cid:54) = k = c ul (cid:54) = k ) = ω u for c k = c uk . The important parameter κ k denotes the impact that a change of concentration c k has on the respective transition rate, i.e. it is a measure15/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 3: Computer simulations of vesicle patterns on real time scale (in min) and spatial scale(scale bar: 5 µ m ). (top) unperturbed cell: { Ω , Ω , Ω F , Ω } = { . , . , . , . } , (bottom)perturbed cell: { Ω , Ω , Ω F , Ω } = { . , . , . , . } . The histograms show the evolution of thevesicle size distribution together with the fitted log-normal distribution (red line). Values for µ and σ are given in Table 5.of sensitivity toward that particular protein. Of course, κ k = κ k ( θ (cid:48) | θ ) in full notation, i.e. thevalue does not only change across proteins, but the impact also changes for different transitions.16/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 Putting eqs. (12), (13) together, we arrive at: ∆ ω k = ω ( c k , c l (cid:54) = k = c ul (cid:54) = k ) − ω u = ω u κ k c uk ( c k − c uk ) (14)In many experimental cases, as e.g. in RNAi screens, the perturbation of a protein concentrationleads to c k → , i.e. we are interested in the dynamics in the absence of a given protein. Withthis assumption we finally have: ∆ ω k ( θ | θ (cid:48) ) = − κ k ( θ | θ (cid:48) ) ω u ( θ | θ (cid:48) ) . (15)This allows us to relate two different dynamical scenarios and their respective outcome in termsof the vesicle patterns: (i) the unperturbed scenario with experimentally known concentrationsand known transition rates, and (ii) the perturbed scenario , where different proteins may beabsent.As an example, let us investigate how changes in the concentration c k of the protein k = CDK8 affect the transition rates ω CDK8 (1 | and ω CDK8 (0 | . Given the parameters in Figure 3, thereduced transition rates return Ω CDK81 − Ω u = 0 . , where Ω CDK8 i refers to the case where theprotein concentration c k = c CDK8 = 0 , whereas Ω ui refer to the unperturbed case. Dividing eqn.(15) by ω b (0 | , we find ∆Ω CDK81 = Ω CDK81 − Ω u = 0 . 375 = − ω u (1 | ω b (0 | · κ CDK8 (1 | − Ω u · κ CDK8 (1 | , (16)from where it follows that κ CDK8 (1 | ≈ − . , which describes the sensitivity toward changesin the concentration of CDK8.Knowing the difference ∆Ω CDK81 , the transition rates ω k are not fully determined because of Ω u = ω u (1 | /ω b (0 | , Ω CDK81 = ω CDK8 (1 | /ω b (0 | . Hence, we can now discuss two differentcases which refer to two hypotheses about the transport of vesicles towards microtubule minus-ends.The first hypothesis of our model states that the transition rate ω CDK8 (1 | in cells that aresilenced for CDK8 reads: ω CDK8 (1 | 0) = ω u (1 | · exp (cid:18) − . · c CDK8 − c u CDK8 c u CDK8 (cid:19) . (17)Assuming c CDK8 ≈ , it follows, that ω CDK8 (1 | is increased by a factor of approximately 1.5with respect to ω u (1 | . This means that silencing the protein CDK8 increases the transition ofvesicles freely diffusing in the cytosol to vesicles being transported towards the MT-plusends. Wethus hypothesize that CDK8 is either directly or indirectly involved in the docking of vesicles onmicrotubule filaments via kinesins. 17/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 On the other hand, κ CDK8 (0 | can be obtained by the following relation: Ω CDK − Ω u = 0 . 375 = ω b (1 | ω CDK (0 | · κ CDK8 (0 | CDK · κ CDK8 (0 | , (18)from where it follows that κ CDK8 (0 | ≈ . . The second hypothesis of our model states thatthe transition rate ω CDK8 (0 | in cells that are silenced for CDK8 reads: ω CDK8 (0 | 1) = ω u (0 | · exp (cid:18) . · c CDK8 − c u CDK8 c u CDK8 (cid:19) . (19)Setting c CDK8 ≈ , we find that ω CDK8 (0 | is decreased by a factor of approximately 0.75 withrespect to ω u (0 | . This means that silencing the protein CDK8 leads to a decrease in transitionof vesicles transported towards the MT-plusends to vesicles freely diffusing vesicles. We thushypothesize that CDK8 may as well be involved in the undocking of vesicles bound to MTfilaments via kinesins.In a similar manner, we could estimate the concentration dependence of other transition rates,using the scaled transition rates (Ω CDK , Ω CDK F , Ω CDK ) and (Ω u , Ω uF , Ω u ) . This enables us tofurther develop a number of hypotheses regarding the effect of silencing the protein CDK8 onthe intracellular transport. Genomic and pharmaceutical research nowadays heavily relies on systematic screens in which theperturbation or silencing of specific proteins affects the abundance of vesicle patterns observedin a population of cells. The study of vesicle pattern formation thus is important to improveour understanding of the function of genes. To learn how vesicle patterns are formed, we haveset up an agent-based model of intracellular transport inside a single cell. Agents representvesicles which move either by diffusion in the cytosol or are transported along the cytoskeletalfilaments through molecular motors. Vesicles further interact with other vesicles by fusion orfission, or they interact with the cytoskeletal filaments, the cell membrane and the nucleus. Thisinteraction is controlled and regulated by specific proteins which are synthesised inside the cell bytranscriptionally active genes. The activity of these genes thus represents the control parametersof our system.Treating vesicles as Brownian agents with an internal degree of freedom allows to formallyderive a model that captures all relevant processes in the formation of vesicle patterns. Five18/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 different values of the internal degree of freedom define the vesicle’s different modes of activity.Transitions between these modes occur as stochastic events related to the binding of proteinsto the vesicle’s coat, or to signalling events. Proteins involved in those processes can controlthe vesicle’s activity which in turn determines the process of pattern formation. We thereforeassume that the transition rates between different value of the internal degree of freedom dependon the concentration of a cell’s proteins. To determine the precise transition rates as functions ofprotein concentrations would require the full knowledge about the regulatory structure of genenetworks. We simplified this situation by assuming that the transition rates can be decomposedinto a product of the unperturbed transition rates and an exponential function depending onthe concentration of single proteins. We further simplified our model by assuming that the cellmembrane, the cytoskeletal filaments and the nucleus are static and that the diffusion coefficientand velocities along cytoskeletal filaments are constant parameters of the model. In contrast,the transition rates represent the free parameters of our model and have been varied. We havereduced the number of 10 original transition rates to 4 scaled transition rates. This introducedan ambiguity in the interpretation of the simulation results with respect to the original transitionrates. We could cope with this situation by deriving two complementary hypothesis about theinfluence of the proteins involved, which could be tested experimentally. In order to relate our computer simulations to reality, experimental data are needed to calibratethe simulations. Whenever available, baseline values for the transition rates obtained from exper-iments have been included. This allowed us to generate patterns on real time and spatial scales.From every simulated pattern, four features can be obtained: size, relative distance to nucleus,number of vesicles within a fixed radius around each vesicle and number of vesicles per cell area.To measure the dissimilarity between the simulated pattern and the experimentally observed one,it is useful to compute the symmetrized Kullback-Leibler divergence of the two correspondingvesicle feature distributions [1], to find out for which parameters the simulated pattern provides aminimum divergence to the experimentally observed vesicle patterns. An interpretation of thesefindings in comparison with hypotheses generated from the computer simulations allows to drawconclusions about the underlying processes, in particular about the role of the genes involved.A comparison of simulated and experimentally measured vesicle patterns also involves a dimen-sional problem: our simulations are performed in 2d, whereas experimental patterns result fromvesicle motion in 3d. In principle, this would require to correct for the distances as consequenceof the projection from 3d to 2d. Since mammalian cells are relatively flat with the exception ofthe nucleus area, we have omitted this correction. But we considered the fact that vesicles couldpass below/above each other by not assuming mutual exclusion in our 2d simulations.19/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 We emphasize again that our computer simulations lead to testable hypotheses about the influ-ence of specific genes, such as CDK8, on the formation of vesicle patterns. One future applicationof our model of intracellular transport concerns the silencing of multiple genes. Let us denotewith i and l two genes which are silenced, then according to ansatz in eqn. (13) we find for thetransition rates ω i ; l ( θ (cid:48) j | θ j ) = ω ( θ (cid:48) j | θ j ) · exp (cid:18) κ i ( θ (cid:48) j | θ j ) c i − c ui c ui + κ l ( θ (cid:48) j | θ j ) c l − c ul c ul (cid:19) . (20)Once we have determined κ i ( θ (cid:48) j | θ j ) and κ l ( θ (cid:48) j | θ j ) from single gene silencing experiments and theirrelated simulations of our model of intracellular transport, we can predict the resulting patternand dynamics of the combined silencing of these two genes. However, (eq. 20) is only valid, ifgenes i and l are not in the same (regulatory) pathway. Such a prediction represents an in-silicoexperiment and can be tested experimentally. Acknowledgment M.B. could benefit from numerous stimulating discussions with Lucas Pelkmans. References [1] Birbaumer, M. (2010). Statistical analysis and modeling of endocytic vesicle pattern forma-tion . Ph.D. thesis, ETH Zurich, Disseration Nr. 19322.[2] Conner, S.; Schmid, S. (2003). Regulated portals of entry into the cell. Nature ,37–47.[3] Dautry-Varsat, A. (1986). Receptor-mediated endocytosis: the intracellular journey of trans-ferrin and its receptor. Biochimie , 375–81.[4] Dinh, A.-T.; Pangarkar, C.; Theofanous, T.; Mitragotri, S. (2006). Theory of Spatial Pat-terns of Intracellular Organelles. Biophysical Journal , L67–L69.[5] Dinh, A.-T.; Pangarkar, C.; Theofanous, T.; Mitragotri, S. (2007). Understanding Intracel-lular Transport Processes Pertinent to Synthetic Gene Delivery via Stochastic Simulationsand Sensitivity Analyses. Biophysical Journal , 831–846.[6] Doherty, G.; McMahon, H. (2009). Mechanisms of Endocytosis. Annual Review of Biochem-istry , 857–902. 20/21 irko Birbaumer, Frank Schweitzer:Agent-Based Modeling of Intracellular Transport European Physical Journal B vol. 82, no. 3-4 (2011) pp. 245-255 [7] Ebeling, W.; Schweitzer, F. (2003). Self-Organisation, Active Brownian Dynamics, andBiological Applications. Nova Acta Leopoldina NF 88(332) , 169–188.[8] Ebeling, W.; Schweitzer, F.; Tilch, B. (1999). Active Brownian particles with energy depotsmodeling animal mobility. Biosystems , 17 – 29.[9] Fiore, P. D.; Camilli, P. D. (2001). Endocytosis and Signaling: An Inseparable Partnership. Cell , 1 – 4.[10] Garcia, V.; Birbaumer, M.; Schweitzer, F. (2011). Testing an agent-based model of bacte-rial cell motility: How nutrient concentration affects speed distribution. European PhysicalJournal B , 235–244.[11] Hannon, G. (2002). RNA interference. Nature , 244–251.[12] Mach, R.; Schweitzer, F. (2007). Modeling Vortex Swarming In Daphnia. Bulletin of Math-ematical Biology , 539–562.[13] Marsh, M.; Helenius, A. (2006). Virus entry: open sesame. Cell , 729–740.[14] Mercer, J.; Helenius, A. (2008). Vaccinia Virus Uses Macropinocytosis and ApoptoticMimicry to Enter Host Cells. Science , 531–535.[15] Roth, M. (2006). Clathrin-mediated endocytosis before fluorescent proteins. Nat Rev MolCell Biol , 63–68.[16] Schweitzer, F. (2003). Brownian Agents and Active Particles: Collective Dynamics in theNatural and Social Sciences . Springer, Berlin.[17] Schweitzer, F.; Garcia, D. (2010). An agent-based model of collective emotions in onlinecommunities. European Physical Journal B77(4)