C.M.'s Science Blog

CM shorts: Ideas between science and philosophy

2012-01-05T11:06:00.002+01:00

I have started an additional blog on Tumblr, called CM shorts, which is devoted to more philosophical topics, such as creativity, free will, collaboration, and so on. Each of these short posts presents one small idea at a time, and the order of the ideas is not very systematic. Yet, I hope the blog as a whole will eventually represent a more or less coherent view point.

Reconstructing fiber networks from confocal image stacks

2011-11-17T07:41:00.001+01:00

We present a numerically efficient method to reconstruct a disordered network of thin biopolymers, such as collagen gels, from three-dimensional (3D) image stacks recorded with a confocal microscope. Our method is based on a template matching algorithm that simultaneously performs a binarization and skeletonization of the network. The size and intensity pattern of the template is automatically adapted to the input data so that the method is scale invariant and generic. Furthermore, the template matching threshold is iteratively optimized to ensure that the final skeletonized network obeys a universal property of voxelized random line networks, namely, solid-phase voxels have most likely three solid-phase neighbors in a $3\times3$ neighborhood. This optimization criterion makes our method free of user-defined parameters and the output exceptionally robust against imaging noise.

Article in ArXiv

Poresizes in random line networks

2011-10-20T11:56:00.001+02:00

Many natural fibrous networks with fiber diameters much smaller than the average poresize can be described as three-dimensional (3D) random line networks. We consider here a `Mikado' model for such systems, consisting of straight line segments of equal length, distributed homogeneously and isotropically in space. First, we derive analytically the probability density distribution $p(r_{no})$ for the `nearest obstacle distance' $r_{no}$ between a randomly chosen test point within the network pores and its closest neighboring point on a line segment. Second, we show that in the limit where the line segments are much longer than the typical pore size, $p(r_{no})$ becomes a Rayleigh distribution. The single parameter $\sigma$ of this Rayleigh distribution represents the most probable nearest obstacle distance and can be expressed in terms of the total line length per unit volume. Finally, we show by numerical simulations that $\sigma$ differs only by a constant factor from the intuitive notion of average `pore size', defined by finding the maximum sphere that fits into each pore and then averaging over the radii of these spheres.

Article in ArXiv

Rendering of Latex-Formula broken

2010-10-13T14:48:00.001+02:00

Sadly, I discovered only recently that the Latex equations in this blog are not displayed any more. God knows how long this situation has already prevailed ! This shows, once again, that one cannot really trust in free services, such as "YourEquation.com". My last PDF backup of the blog was, unfortunately, in July 2009, but this is better than nothing (-: I am now thinking to write future blogs as PDF documents.

Poster on "Tumor Cell Migration"

2010-04-03T15:07:00.000+02:00

We have presented a poster at the Spring Meeting of the German Physical Society (DFG) in Regensburg.

On the apparent correlations between diffusivity D and power-law exponent b

2010-03-12T10:53:00.021+01:00

When analyzing the trajectories of cytoskeleton-bound beads or whole cells, one frequently finds MSD curves as a function of lag time that can be fitted to a power-law within a broad temporal range.

To demonstrate such a power-law regime, one plots the MSD double-logarithmically. To show that analytically, one defines an arbitrary unit of time t0 (for example, t0 = 1 min) and length r0 (for example, r0 = 1 um), in order to make the lagtime and the displacement dimensionless. Then, the MSD curves can be locally written as


\left( \Delta R / r_0 \right)^2 = D \cdot \left( \Delta t / t_0 \right) ^b

The two dimensionless parameters D and b are called the (apparent) "diffusivity" and the "power-law exponent", respectively. The logarithm of the above equation reads


\log\left( ( \Delta R / r_0 )^2 \right) = \log(D) + b \cdot \log( \Delta t / t_0)

Defining


y = \log\left( ( \Delta R / r_0 )^2\right)

and


a =  \log(D)

and


x =  \log( \Delta t / t_0)

we obtain a linear relation for our logarithmic variables:


y = a + b x.

Typical experiments yield a whole bunch of MSD-curves, corresponding to a set of N parameter pairs {(a_i,b_i)}. The a_i and b_i are fluctuating, with the a_i often being normally distributed. When the value-pairs (a_i,b_i) are plotted as a point cloud in the a-b plane, one frequently finds correlations, such as high a-values (=logarithmic diffusivities) coming together with small b-values (=power-law exponents).

However, note that these statistical properties depend on the choice of the length and time units that have been chosen to make the variables dimensionless. In the logarithmic domain, the a-value is just the y-intercept of the linear curve, a=y(x=0). In the original domain, this implies


D = e^a = e^{y(x=0)} =\Delta R^2(\Delta t=t_0) / r_0^2.

Therefore, the distribution P(D) will depend on the parameters t_0 and r_0. For an extreme example, imagine a double-log plot with a bunch of straight MSD lines that all intersect at some point. If we evaluate the distribution P(D) precisely at the lagtime of the crossing point, we will obtain a delta-distribution ! Of course, in reality the lines will not all cross at one point, yet they might approach each other closely within some finite spot.

It is convenient to further analyze the situation in the logarithmic domain, and later transform back to the original domain.

The problem is: Given a bunch of N straight lines,


y_i = a_i + b_i \cdot x,

which x=x_opt would be best suited to evaluate the distribution of the y_i(x=x_opt), or later the correponding D_i = e^{y_i(x=x_opt)} ?

I suggest to choose the x_opt where the variance of the y_i becomes minimal:


var\{y_i(x=x_{opt})\} = min.

A straight-forward calculation shows that this x_opt is given by


x_{opt} = - cov\{a_i , b_i\} / var\{b_i\},

where


cov\{a_i , b_i\} = < (a_i-\overline{a}) (b_i-\overline{b}) >_i

and


var\{b_i\} = < (b_i-\overline{b})^2 >_i

with the notation


\overline{c} = (1/N)\sum_{i=1}^N c_i = < c_i >_i

To demonstrate the effect of x_opt on the statistical properties of D and b, I have generated some artificial set of MSD-curves (Actually it consists of N=100 curves, but only 10 are shown for clarity. The slopes/power-law-exponents are unrealistically large, never mind. The units of time and length were 1 min and 1 um, respectively):

If we evaluate the statistics of the (D_i,b_i) at lagtimes dt=0.1min and dt=100min, we obtain in the D-b-plane point clouds shown below in green and red colors, respectively. When using instead the optimum x_opt=0.339, corresponding to t_opt=t0*e^{x_opt}=2.182, the blue point cloud is obtained.

Obviously, most of the apparent correlations have disappeared (Well, not quite: A closer inspection with higher N shows that even at t_opt the variance of the D_i changes systematically with b, and vice versa).

High kurtosis by sum-of-exponential distributions

2009-09-21T16:56:00.006+02:00

The exponential probability distribution


P(x) \propto e^{-\lambda x}

has an excess kurtosis of 6 for arbitrary decay constants.

Different values of the kurtosis can be obtained for distributions that are superpositions of exponential distributions with different decay constants.

As a test, the distribution


P(x) \;\propto\; 10 e^{- 10 x}\; + \; e^{-x}

has been sampled numerically. The resulting set of random numbers has then been analyzed statistically:

The numerical kurtosis of the superposition is 17.

This could be relevant for the "spider web model of cytoskeletal fluctuations": An elementary remodeling step in a remote fiber causes a smaller shift of the bead, i.e. the corresponding probability distribution of shifts has a steeper decay. The total shift distribution is a superposition of the different fiber contributions.

[10] Cell Invasion: Summary of preliminary results

2009-09-15T15:05:00.015+02:00

A. Modelling the cell invasion process

We have developed a set of mathematical models for the invasion process of tumor cell populations into a half-space of collagen gel. In particular, we aimed for a quantitative understanding of the characteristic shape of the invasion profiles, their temporal evolution, and their dependence on the initial cell surface density.

On the microscopic scale, collagen is a highly inhomogeneous fiber network. The detailed procedure by which individual tumor cells migrate through this porous network (involving steps such as finding adhesion ligands, forming and disintegrating focal adhesion contacts, up- and down-regulating acto-myosin traction forces, ..) is not well understood at present. Moreover, it depends on experimentally inaccessible local conditions. It is therefore reasonable to describe cell migration as a stochastic process, i.e. essentially as a random walk. The effective diffusion constant of this random walk summarizes the complex interactions of the cell with its surrounding material in a coarse-grained way. For simplicity, we have further assumed that the collagen gel is statistically homogeneous and isotropic, i.e. the diffusion of (isolated, non-cooperating) cells of a certain type is equally fast at any position within the gel and does not depend on the spatial directions.

B. Figures of merit for testing invasion models

The most detailed description of a (real or simulated) invasion experiment consists of a complete 3D trajectory


\vec{R}_n(t)=(x_n(t),y_n(t),z_n(t))

for each individual cell n. So far, we have mainly focused on the z-coordinates.

In analogy to the time-averaged "mean squared displacement", a standard tool for analyzing stationary random walks, we can define a population-averaged "mean squared invasion depth":


\left\langle z^2 \right\rangle (t) = (1/N) \; \sum_{n=1}^N (z_n(t))^2,

where n runs over all N cells of the population.

The time-dependent z-density distribution (unit: 1/m)) is defined by


n(z,t) = \sum_{n=1}^N \delta(z-z_n(t)).

Since n(z,t) is usually a very noisy quantity, we work in the following with the (equivalent) cumulative probability distribution


P(z,t) = \int_z^\infty n(z,t) dz\; /\; \int_0^\infty n(z,t) dz.

This dimensionless quantity can be interpreted as the probability to find a cell at depth z, or deeper. For any fixed time t_0, the function P(z,t_0) starts with the value 1 at the top surface (z=0) of the collagen slice and monotonically decreases towards zero for larger depths z.

C. Experimental key features

Measured cumulative probability distributions reveal a large variability, depending on the cell type and on the material parameters of the used collagen gel. However, there are a few general properties that we observed in almost all cases investigated so far:

1.) A typical distribution, for fixed time t_0, consists of up to three distinctive layers, or zones. These zones are most easily visable in a semi-logarithmic plot of P(z,t_0).

2.) Within a narrow layer close to the surface, P(z,t_0) is rapidly decaying, i.e. according to a faster-than-exponential law.

3.) Within a broad intermediate zone, the cummulative probability is decaying exponentially, P(z,t_0) ~ e^(-z/z_0) with a characteristic length scale z_0.

4.) For very large depths, at the "front zone", P(z,t_0) decays apparently according to a faster-than-exponential law. Note that, for all finite populations, there is always a single cell at the foremost front of the distribution. Statistics is becoming extremely poor close to that region.

6.) As time passes, the characteristic length scale z_0 of the exponential zone increases.

7.) The mean squared invasion depth of the cell population is increasing with time t. The functional dependence of on t is non-linear, indicating an anomaleous diffusion process.

8.) The invasion profiles P(z,t) depend strongly on the 2D density of tumor cells initially plated onto the gel surface:

9.) For very small population densities, the surface layer, marked by the rapid drop of cumulative invasion propability, disappears.

10.) For increasingly higher densities, the drop of cumulative probability at small depths becomes more pronounced, indicating a dense aggregate of tumor cells immediately below the surface. Beyond this aggregate, only a relatively small fraction of cells is invading into the deep bulk of the gel. However, the characteristic length scale z_0 of the invaded fraction tends to increase with initial cell density.

D. Various investigated models

We have considered a variety of different invasion models and explored their ability to account for the above experimental key features.

The most basic possibility would be a homogeneous population of non-cooperative tumor cells, where all individuals migrate with the same effective diffusion constant, independently from each other. However, the simple Gaussian invasion profiles, resulting from standard diffusion into a half-space, are clearly inconsistent with the complex layered structure of the measured P(z,t) distributions, in particular with the zone of exponential decay (experimental feature 3).

Heterogeneous cell populations with a spectrum of vastly different diffusion constants could produce close-to-exponential cumulative probabilities. Yet, this requires fine-tuning of the assumed parameter distribution. In addition, heterogeneity cannot account for the observed cell density dependence.

It is well-known that standard diffusion combined with a process of reverse drift, back to the surface, at constant velocity, produces exponential profiles. One then has to explain the biological cause of this slight trend in the random walk of the tumor cells. Among the possibilities are mesoscopic spatial gradients in the distribution of local properties (topological, rheological, or chemical) within the gel. Again, such external causes cannot account for the observed cell density dependence.

We have also considered a model in which the cells themselves emit a chemo-attractant, which independently diffuses within the gel and degrades at a certain rate. For certain parameter ranges, the resulting invasion profiles are indeed consistent with the observations. However, the values of the diffusion constants required for a good fit to the data turned out to be not realistic.

E. The "Density-Dependend Diffusion" model

At present, our best model for the invasion process is based on the assumption of a situation-dependent diffusion constant of the tumor cells. This idea is motivated by the observation of the cell aggregate (experimental feature 10) that formes in a narrow layer below the gel surface when the initial cell surface density was large. Within this layer, neighboring cells are very close to each other. The aggregate remains rather stable over long times, meaning that the diffusion constant of cells in an aggregated state is small. On the other hand, once individual cells escape from the aggregate into the bulk of the gel, their diffusion constant becomes much higher, indicated by the large characteristic length scale of the exponential zone.

We therefore tested a simple model in which the individual cells switch their diffusion constant between a small and a large value, depending on the presence or absence of at least one neighboring cell within a certain detection range. The full model has been implemented as a multi-dimenional Monte-Carlo simulation. Additionally, in a mean field approximation, the corresponding 1D Fokker-Planck equation has been solved numerically.

The three profiles shown in the figure below correspond to three separate invasion experiments, differing only in the initial cell surface density (The red, green and blue curve correspond to a total number of 33, 498 and 4348 cells in the total field of view). In all three cases, the tumor cell positions within the gel have been measured after the same time delay (relative to planting the cells onto the gel surface). The characteristic density dependent changes of the invasion behavior, especially the features 9 and 10 listed above, are clearly visible.

The next figure shows the results of corresponding computer experiments, obtained by the Monte Carlo simulation of the DDD model. No attempt has been made to fit the parameters to the specific experiments. Nevertheless, all major features and trends are qualitatively reproduced:

[9] Cell-Invasion: Implementation of the DDD model

2009-08-12T13:53:00.012+02:00

It is possible to solve the Fokker-Planck (FP) equation of post [8] numerically by direct temporal forward integration, as done before.

Let us discretize space as


z_i = i \; \Delta z

and time as


t_n = n \; \Delta t.

Let \mu be the fraction of density transfered from site i to site i+1 during one time step, so that the corresponding current would be


J_{i,i\!+\!1} = \mu \; \rho_i

and let us denote the density in site i at time n as


\rho^{n}_i = \rho(z_i,t_n).

Except for the most left (i=0) and the most right (i=i_mx) sites (where boundary conditions have to applied properly), the discretized diffusion equation reads


\rho^{n+1}_i  = (1-2\mu) \rho^{n}_i + \mu (\rho^{n}_{i-1} + \rho^{n}_{i+1}).

Note that the physical diffusion constant is related to the discretization parameters as


D = \frac{\mu \Delta z^2}{\Delta t}.

To avoid numerical instabilities, one can set \mu=0.1 and choose a proper time step.

We shall later apply the above method to the DDD model. In this case, \mu is locally replaced by


\mu \rightarrow \mu_i = \mu(\rho^{n}_i).

As a simple test of the method, we start with ordinary diffusion in a box of width 100 with hard walls at z=0 and z=100. The initial density at t=0 is a narrow Gaussian, centered in the middle of the box. The diffusion constant was set to 10.

In a semilog plot, one gets the expected evolution of density:

The variance of z-distributions follows initially the behaviour of free diffusion. Later it saturates because of the confinement in the box:

Next we turn to the model with density dependent diffusion constant (DDD model). The initial density is now localized in the z=0 layer. The evolution of the particle densities, as resulting from numerical integration of the FP equation, looks as follows:

In the semilog. plot one can see the characteristic bi-phasic behaviour and the almost exponential regime at intermediate z.

Here the corresponding cummulative probabilities (= invasion profiles):

The variance as a function of time ...

... starts sub-diffusively and later appears to become slightly super-diffusive.

[8] Cell-Invasion : Clustering : Density Dependent Diffusion (DDD) Model

2009-08-12T13:08:00.007+02:00

It would be convenient to formulate the clustering effect in the form of a 1D Fokker-Planck equation that can be quickly solved numerically. For that purpose, remember that the main point of the clustering effect is the simultaneous presence of fast and slow diffusing particles at any depth z. Assume we know the fraction P_NC of particles at z which have no closeby partner within their interaction range r_c. Then we can define a local effective diffusion constant as


D_{eff} =  P_{NC} D_{fast} + (1-P_{NC}) D_{slow}.

Consider a small stripe of material around depth z. Let the local 3D particle density within this z-stripe be rho. We neglect any possible correlations between the positions of different particles within the stripe, i.e. we assume spatial Poisson statistics with average density rho. Then the fraction P_NC is the probability that, for an arbitrary particle, the nearest neighbor is found in a distance larger than r_c. The average number of particles in a sphere of radius r_c is given by


\overline{n} = \rho \;\frac{4}{3}\pi r^3_c.

According to Poisson statistics, the probability that n=0 is


P_{NC} = e^{-\overline{n}} = e^{ - \rho \;\frac{4}{3}\pi r^3_c} = e^{- \rho / \rho_c},

with


\rho_c = \frac{3}{4 \pi}\; r^{-3}_c.

Thus we can define a density dependent diffusion constant:


D_{eff}(\rho) \;= \; (e^{-\rho / \rho_c})  D_{fast} \;+\; (1-e^{-\rho / \rho_c}) D_{slow}.

The figure below shows an example:

Using this, we can write our model as a non-linear Fokker-Planck equation for the position and time-dependent particle density:


\frac{d}{dt} \rho(z,t)\; = \;D_{eff}(\rho(z))\; \frac{\partial^2}{\partial x^2} \;\rho(z,t).

[7] Cell-Invasion : Clustering : Bi-phasic invasion profiles

2009-08-10T17:24:00.008+02:00

Another interesting feature of the clustering effect (compare posts [4] and [5]) is the emergence of two relatively distinct phases in the invasion profiles:

Close to the cell monolayer at z=0, the cell density is high, so that clusters have a high probability. This leads to a small average diffusion constant and, thus, a weakly invasive phase close to the surface.

At some point, the absolute cell density and the corresponding cluster probability fall below the critical value and most cells switch to a high diffusion constant. A strongly invasive phase is formed further in the bulk.

The following MC simulation demonstrates this effect:

Remarkably, many of the measured invasion profiles show qualitatively such a bi-phasic behaviour:

Note that for technical reasons the coordinate of the monolayer is not at z=0. No attempt has been made to fit the data.

--------------

Some side remark: As mentioned already in [6], the value P_cum[z=1] can be interpreted as the invaded fraction of cells. When we zoom into the simulation result above, we can see that this fraction increases from about 0.63 to 0.86 within the time interval between t=20 and t=100 units.

[6] Cell-Invasion : Slow Layer Effect

2009-08-10T13:06:00.021+02:00

The cell clustering effect, as described in posts [4] and [5], is not easy to treat analytically. A much simpler, but somewhat similar model is the "slow layer effect":

The particles (tumor cells) undergo standard diffusion (with rates Rb) in the whole bulk of the material. Only in the z=0-layer the diffusion constant (i.e. the transfer rates for intra layer motion and for steps out of the layer) is reduced to a value Rl smaller than Rb.

Before we go into simulations, a simple one-dimensional consideration shows that the abrupt change of rates at z=0 implies a sudden drop of particle density. We have for the temporal change of density \rho_0 at z=0:


\frac{d}{dt}\rho[0] = -\rho[0] R_l + \rho[1] R_b

If we assume a quasi-steady state, where rho[0] and rho[1] change very slowly with time, we obtain


\frac{\rho[1]}{\rho[0]} = \frac{R_l}{R_b},

and thus rho[1] is much smaller than rho[0], if R_l is much smaller than R_b.

The following figure shows simulation results for the invasion profiles (cummulative probabilities) in a linear plot. The ratio R_l/R_b was 1/1000. The curves correspond to different times (100,200,...,500 units):

One can clearly see the drop of density immediately above the z=0-layer. Note that P_cum[1] acan be interpreted as the fraction of particles that have separated from the boundary layer and invaded the bulk. We might call this quantity the "invaded fraction". It is about 0.4 in the magenta-colored profile above.

A semilog plot of the same results shows that effects of the slow layer range far into the bulk:

The shape of P_cum(z) is slightly different from the case of standard diffusion in [4]. The initial "slope" of the curves appears steeper.

Actually, if one removes the single data point at z=0, P_cum(z) can very well be fitted (in the whole bulk z=1...zMax) by shifted Gaussians, using 3 free parameters P_0, z_0 and \sigma:


P_{cum}(z>1) = P_0 \; e^{-(z+z_0)^2/\sigma^2}

Basically this means that the curves are approximately factors of exponentials and Gaussians. For not too large z, the exponential factor will dominate. The exponential range should grow with sigma. Therefore, the larger the bulk diffusion constant, the invasion profiles should appear more and more exponential-like. Hence, even the simple "slow layer" model can to some extent account for exponential invasion profiles.

Is it possible to treat this model analytically in 1D ? I guess that in the limit of smaller and smaller width of the z-slices, the effect of the slow layer translates into a simple boundary value condition for the derivative of rho(z). Except this unusual boundary condition, ones has only to solve a standard diffusion Fokker-Planck equation in the bulk.

I tried a very simple numeric solution of the FP equation, using the methods of an old post (Program: Invade1):

Note that this is the real particle density P(z), not the cummulative invasion profile. The density jump corresponds to the expected value Rl/Rb. One can also see that the initial slope of the bulk density curves is negative.

Finally, the corresponding cumulative invasion profile:

[5] Cell-Invasion : Clustering causes fractional diffusion behaviour ?

2009-08-09T09:22:00.022+02:00

One of the most used quantities for the characterization of random walks is the Mean Squared Displacement (MSD) as a function of lagtime, averaged over absolute time. In the case of cell invasion (which might be non-stationary and therefore not suitable for being analyzed with standard MSD), we can define similar quantities of interest as population averages of (powers of) the z-coordinate versus absolute time. For example, we might consider

the Mean Invasion Depth (MID):


\overline{z}(t)=\left\langle z \right\rangle_{pop}=\frac{1}{N}\sum_{n=1}^N z_n(t),

where index n runs over all N agents,

the Mean Squared Invasion Depth (MSID):


\overline{z^2}(t)=\left\langle z^2 \right\rangle_{pop}=\frac{1}{N}\sum_{n=1}^N z^2_n(t),

the Variance (var):


\mbox{var}(t)=\left\langle (z-\overline{z})^2 \right\rangle_{pop}=\frac{1}{N}\sum_{n=1}^N (z_n-\overline{z})^2(t),

and the Invasion Front Position (IFP):


z_{front}(t)=\mbox{max}\left\{ z_n(t) \right\}_{n=1\ldots N}.

Some of these quantities are evaluated in the following for the cell invasion model presented in the last post [4]. We first look at the case of independent random walkers:

In the double-logarithmic plot, the MSID is found to grow linearly with time, as would be expected for normal diffusion. The MID, accordingly, grows with the square root of time. The front particle is always ahead of the mean, but is advancing with the same square-root law.

Now we turn to the invasion with clustering:

Surprisingly, the MSID grows with time as a fractional, super-diffusive power law ! The MID also grows faster than with the square root of time.

Isn't that counter-intuitive ? How can the occasional slowing down of clustered agents lead, effectively, to a "faster-than-diffusive" invasion dynamics ?

A closer look at the figures reveals that, within the first 100 time units, the normally diffusing agents have invaded - on average - further than the clustering agents, in accordance with common sense. Nevertheless, due to the higher power-law exponent of the clustering agents, there would be a cross-over at some later time. However, the MID and the MSID in the model with clustering seem to show a change of behaviour after about 100 time units. They then seem to return to normal diffusive behaviour, thus avoiding a cross-over.

So the super-diffusion is only a transient effect. It remains surprising, at least for me.

[4] Cell-Invasion : The clustering effect

2009-08-08T19:07:00.010+02:00

Last week I had a nice collaboration with Sebastian Probst. He had a very interesting idea:

Assume the agents have a tendency to form and stabilize clusters: Each time they happen to find another closeby agent on their random walk, they drastically reduce their diffusion constant. They won't stop moving entirely, but let's assume they slow down whenever they are in happy company. How would that clustering effect change the invasion process ?

In order to model this idea, we first start again with the simple diffusion simulation from the last post [3]. In regular time intervalls of 100 units, we count how many agents are found within each z-stripe, resulting in a histogram N(z). From this histogram we compute a cumulative probability distribution P_cum(z), which gives the probability that an agent is found in depth z of the material or deeper. In the following, we call this cumulative distributions the invasion profiles.

One finds for the simple diffusion P_cum(z)-curves which have a negative curvature and approximately Gaussian shapes with a variance that grows with time. Note that this is a semi-log plot, so true Gaussians would show up as downward parabolas.

Now we leave everything as before, except that the clustering effect is additionally implemented. For this purpose, we count the number of direct neighbors of each agent. If this number is larger than zero for certain agents, we reduce the hopping rates of these specific agents by a factor of 10 during the next MC cycle. Again we plot the invasion profiles:

The result is remarkable:

As a consequence of the clustering effect, the profiles are changed from about Gaussian to almost exponential shape (showing up as linear lines in the semi-log plot).

This offers a completely new, simple, and testable mechanism for explaining the mostly exponential invasion profiles observed in the experiments. It does not require any explicit assumptions about the diffusion of chemical signalling molecules. The only parameters are the maximum diffusion constant of the agents, their detection range for sensing closeby fellows, and the reduced diffusion constant in the clustering mode.

Future tests:

The initial density of the cells in the monolayer should drastically affect the invasion process. A small density will reduce the probability of cluster formation and thus drive the system back to the normal, Gaussian diffusion of the agents.

[3] Cell-Invasion : 2D Monte-Carlo simulation

2009-08-08T18:53:00.009+02:00

We have implemented a Monte-Carlo simulation program for studying tumor cell invasion into a half space of complex bio-material (as a part of a new project described in a recent post and its followups). In the following, the tumor cells are treated as active agents in a multi agent system.

Space is discretized as a 2D rectangular grid of sites (x,z). The positive z-direction corresponds to the invasion direction, pointing into the material.

Each site can be occupied by at most one agent. Multiple occupation is forbidden (in a way that is analogous to the Hubbard model used in the theory of metal insulator transitions).

Each site is assigned a potential energy U(z,x), reflecting the difficulty for an agent to occupy this site. It would be possible, in principle, to investigate potential landscapes of arbitrary complexity. Simple choices would be a completely flat potential, or a binary valued potential without spatial correlations, i.e. binary white noise. The two possible values of the local potential could be zero (no obstacle) and infinity (impenetrable obstacle). And so on.

There is a certain probability per unit time that an agent moves from its present site to a neighboring site (This leads formally to a kind of "hopping rate", although the agents are thought to move continuously in reality. It is the discrete, "pixelized" detection of their position that implies discrete hopping steps). The maximum hopping rate is R0. However, when the agent moves from a site of low to a site of high potential energy, the hopping rate is reduced by a Boltzmann factor with some effective temperature (Of course, tumor cell invasion is not driven by thermal fluctuations, but by active motor forces. Nevertheless, the cell migration rate will slow down when the cell tries to move into obstacles, i.e. when it enters high energy sites in the potential landscape. This effect can be conveniently modelled by the Boltzmann factor.)

For the horizontal (z) and vertical (x) edges of the simulation area, boundary conditions can be chosen between hard walls and periodic in our implementation.

Many agents can be in the system simultaneously. They can interact with each other (in order to describe processes such as binding of tumor cells etc.). They can also change the potential landscape (describing, for instance, stigmergic effects such as the protolysis of collagen by chemicals emitted by the tumor cells).

The Monte-Carlo algorithm is "exact": In each cycle, a complete list of all possible processes (all possible hops of all agents, etc.) is generated and the rates of these processes (depending on the current system state) are computed, as well as the sum of all rates, which determines the average waiting time dtMean until the next process happens. The MC algoritm then chooses a random, exponentially distributed waiting time according to dtMean. It also chooses one of the processes from the list, with the correct relative probabilities, according to the rates. The process chosen is executed, the system state is updated, the relevant variables are recorded, and the next cycle starts.

As a simple example, let us consider a MC simulation on a grid of 1000(x) times 100(z) sites. In the x-direction, we chose periodic boundaries, but in the z-direction the boundaries are hard (This is realized by placing an almost infinitely high potential wall just outside and along the west and east edge of the potential landscape). Initially, 1000 cells are placed onto the west surface (z=0), forming a complete monolayer. The cells then start to diffuse randomly (This is simply achieved by a flat potential landscape U(z,x)=const. inside the whole simulation area. Then the hopping rates are everywhere homogeneous and isotropic). The picture below shows a snapshot of the distribution of agents over the simulation area after 100 time units (black circles) and after 500 time units (red squares). One can see how the agents are invading the half space.

[2] Interpretation of the directional correlation function

2009-07-21T14:13:00.016+02:00

In post [1] on random walks with directional persistence, we have defined the correlation function of the direction vectors as

C_{mn}=C_{m-n}=\left\langle \vec{e}_m  \vec{e}_n \right\rangle = \left\langle \cos(\phi_m-\phi_n) \right\rangle,

where the brackets denote the ensemble average. We note that

 C_1 = \left\langle \cos(\phi_{m+1}-\phi_m) \right\rangle.

The quantity

 \Delta \Phi_1=\phi_{m+1}-\phi_m

is the so-called turning angle, i.e. the angle between the directions of two successive moves of the random walker. Therefore, the value of correlation function at lagtime 1 is just the average of the cosine of the turning angle. If the probability distribution of the turning angle is denoted by

 P(\Delta \Phi_1),

we thus have

 C_1 = \int_{-\pi}^{\pi} d\Delta \Phi_1 \;P(\Delta \Phi_1) \; \cos(\Delta \Phi_1).

In earlier publication s, researchers have defined an index of persistence by

 p_{\Phi_1} = 2 \left( \int_{-\pi/2}^{\pi/2} d\Delta \Phi_1 \;P(\Delta \Phi_1)\right) -1.

This can be rewritten in the form

  p_{\Phi_1} = \int_{-\pi}^{\pi} d\Delta \Phi_1 \;P(\Delta \Phi_1)\; g(\Delta \Phi_1),

where the weight function g() has the constant value +1 in the central interval [-pi/2...+pi/2] and the value -1 in the remaining intervals [-pi...-pi/2] and [+pi/2..+pi]. It is clear that this weight function can be equally well be replaced by the smooth cosine function above. Therefore, we learn that the index of persistence is basically the value of the directional correlation function at lagtime 1:

 p_{\Phi_1} \approx C_1.

In an analogous way, the correlation function at higher lagtimes N>1 is just the index of persistence at the corresponding lagtime N:

 p_{\Phi_N} \approx C_N,

provided that the turning angles are properly defined as

 \Delta \Phi_N = \phi_{m+N}-\phi_m.

[1] MSD of random walks with given SWD and directional correlations

2009-07-20T15:14:00.071+02:00

We consider a stationary random walk of a particle in 2D. The position of the particle is measured at regular time intervalls

t_n=n\Delta t \; \mbox{with} \; n=0,1,\ldots

After each intervall, the particle has moved by one step, characterized by a positive step width and a direction (unit) vector:

\Delta\vec{R}_n = s_n\;\vec{e}_n = s_n \; \left( \cos(\phi_n) , \sin(\phi_n) \right) .

Let the probability distribution function of the step widths (abbreviated by SWD) be given by

P(s) = \mbox{Prob}(s_n=s)

with finite mean and variance and assume that the successive step widths are statistically independent random variables.

\left\langle (s_m-\overline{s})(s_n-\overline{s})\right\rangle = \overline{s^2}\delta_{mn}.

Let the correlation function of the direction vectors be defined as

C_{mn}=\left\langle \vec{e}_m  \vec{e}_n \right\rangle = \left\langle \cos(\phi_m-\phi_n) \right\rangle,

where the brackets denote the ensemble average. Assume that the width and direction of each step are statistically uncorrelated.

\left\langle (s_m-\overline{s})\vec{e}_n\right\rangle = \vec{0}.

The position of the particle at time step N is given by

\vec{R}_N = \sum_n \Delta\vec{R}_n.

What is the ensemble averaged mean squared displacement (MSD) of the particle ?

We have

\left\langle \vec{R}^2_N  \right\rangle = \left\langle  \sum_m \sum_n \Delta\vec{R}_m \Delta\vec{R}_n \right\rangle

=\sum_{m,n} \left\langle s_m\vec{e}_m\;s_n\vec{e}_n \right\rangle

Due to the independence of step lengths and directions this factorizes to

=\sum_{m,n}  \left\langle s_m s_n \right\rangle \; \left\langle \vec{e}_m \vec{e}_n  \right\rangle

=\sum_{m,n} \left\langle s_m s_n \right\rangle \; C_{mn}.

We next write the step widths in the form of average plus fluctuation:

s_m = \overline{s} + \Delta s_m,

understanding that the fluctuations are zero mean:

\left\langle \Delta s_m \right\rangle = 0.

This yields

\left\langle \vec{R}^2_N  \right\rangle = \sum_{m,n} C_{mn} \left\langle (\overline{s} + \Delta s_m) (\overline{s} + \Delta s_n) \right\rangle.

Under our statistical assumptions, we have

\left\langle (\overline{s} + \Delta s_m) (\overline{s} + \Delta s_n) \right\rangle = \overline{s}^2 + \left\langle \Delta s_m \Delta s_n \right\rangle =   \overline{s}^2 + \overline{s^2}\;\delta_{mn},

Note that here

 \overline{s^2}

stands for the variance of the fluctuation of the step widths. We obtain

\left\langle \vec{R}^2_N  \right\rangle = \overline{s}^2\sum_{mn}C_{mn}+\overline{s^2}\sum_m C_{mm}.

Since

 C_{mm}=\left\langle \vec{e}^2_m \right\rangle = 1,

we obtain

 \left\langle \vec{R}^2_N  \right\rangle = \overline{s}^2\sum_{mn}C_{mn}+N\;\overline{s^2}.

Because our random walk was assumed to be a stationary random process, the correlation functions can only depend on time differences,

 C_{mn}=C_{m-n}

, and thus we have

 \left\langle \vec{R}^2_N  \right\rangle = N\;\overline{s^2} + \overline{s}^2\sum_{mn}C_{m-n}.

It is easy to show that

 \sum_{m=1}^N \sum_{n=1}^N f_{(m-n)}=\sum_{k=-N+1}^{N-1}(N-k)  f_k.

Using this, we arrive at

 \left\langle \vec{R}^2_N  \right\rangle = \overline{s^2}\;N + \overline{s}^2\sum_{k=-N+1}^{N-1}(N-k)C_k

 =\overline{s^2}\;N + N\overline{s}^2\sum_{k=-N+1}^{N-1}C_k - \overline{s}^2\sum_{k=-N+1}^{N-1}k C_k.

Since the correlation function is symmetric with respect to k, while k is asymmetric, the last term vanishes, yielding

 \left\langle \vec{R}^2_N  \right\rangle = \overline{s^2}\;N + \overline{s}^2\; N\;\sum_{k=-N+1}^{N-1}C_k.

In the sum, we can further extract the k=0 term

C_0=1

and make use of the symmetry

C_{-k}=C_k

to obtain

\left\langle \vec{R}^2_N  \right\rangle = \overline{s^2}\;N + \overline{s}^2\;N + \overline{s}^2\; 2N\;\sum_{k=1}^{N-1}C_k.

Thus we finally arrive at

\left\langle \vec{R}^2_N  \right\rangle = \left( \overline{s^2}+ \overline{s}^2+ 2\overline{s}^2\sum_{k=1}^{N-1}C_k\right) N.

We see that the MSD of our particle is completely determined by the mean and variance of the step width distribution and by the directional correlation function. All step width distributions with the same mean(s) and var(s) produce the same MSD. Moreover, interesting (non-diffusive) MSD-versus-lagtime relations can only be generated via the directional correlations. We have excluded the possibility of Levi flights and the like by demanding that the SWD is not heavy-tailed.

Note to self: There are at least 4 possibilities how fractional powerlaw MSDs can arise: (1) heavy-tailed step width, (2) heavy-tailed waiting times, (3) persistent directions, or (4) long-time correlated statistical dependencies between the widths and directions of the steps. It might be interesting to follow possibility (4), if I should ever find the time.

We can check a few limiting cases:

(a) Brownian diffusion

In this case

 C_k=\delta_{k0},

so that

 \sum_{k=1}^{N-1}C_k = 0

and, hence,

 \left\langle \vec{R}^2_N  \right\rangle = (\overline{s}^2 + \overline{s^2})\;N,

as it should be for a diffusive process. Note that for an exponential step width distribution,

\overline{s}^2 = \overline{s^2},

so that

 \left\langle \vec{R}^2_N  \right\rangle = 2 \overline{s}^2 N.

(b) Ballistic motion

In this case

 \phi_m=\phi_n=\phi,

so that we get

 \cos(\phi_m-\phi_n)=\cos(0)=1

and so

 C_{mn}=1=C_k

and so

 \sum_{1}^{N-1}C_k = (N-1)

and so

\left\langle \vec{R}^2_N  \right\rangle = \left( \overline{s^2}+ \overline{s}^2+ 2\overline{s}^2 (N-1)\right) N.

In the limit of large N we obtain

\left\langle \vec{R}^2_N  \right\rangle = 2\overline{s}^2 N^2,

i.e. the MSD grows quadratically with time, as expected for ballistic motion.

(c) Exponentially decaying directional correlations

In this case,

C_k=q^{-k}.

For a measuring time intervall dt and a physical decorrelation time \tau, we would have to set

 q=e^{\Delta t/\tau}.

The sum can be written as

\sum_{k=1}^{N-1}C_k = -1 + \sum_{k=0}^{N-1} (q^{-1})^k = -1 +\frac{1-q^{-N}}{1-q},

so that

\left\langle \vec{R}^2_N  \right\rangle = \left( \overline{s^2}+ \overline{s}^2+ 2\overline{s}^2 (-1 +\frac{1-q^{-N}}{1-q}) \right) N.

=\left( \overline{s^2}- \overline{s}^2+ 2\overline{s}^2 \frac{1-q^{-N}}{1-q} \right) N.

Note that for

N\Delta t/\tau<<1

we have

 1-q^{-N}=1-e^{-N\Delta t/\tau}\approx 1-1+N\Delta t/\tau=N\Delta t/\tau=N \ln(q).

So for lagtimes smaller than the decorrelation time, the MSD is ballistic, for larger lagtimes diffusive, as it should be.

(c) Long-time correlated directional correlations

We next assume that for large lagtimes, the correlations decac according to a power law

C_k\rightarrow c_0 N^{-\gamma} = c_0 (t/\Delta t)^{-\gamma}.

To treat this case, we change from discrete time steps to a continuous time by replacing

 N\rightarrow t/\Delta t

and

\sum_{k=1}^{N-1}C_k \rightarrow \frac{1}{\Delta t}\int_{t0}^t C(t^{\prime})dt^{\prime},

to obtain

\left\langle \vec{R}^2(t)  \right\rangle = \left( \overline{s^2}+ \overline{s}^2+ \frac{2\overline{s}^2}{\Delta t}\int_{t0}^t C(t^{\prime})dt^{\prime}\right)t/\Delta t .

In the limit of large lagtimes we have

\int_{t0}^t C(t^{\prime})dt^{\prime} = \int_{t0}^t c_0(t^{\prime}/\Delta t)^{-\gamma}dt^{\prime} \rightarrow (t/\Delta t)^{1-\gamma},

so that

\left\langle \vec{R}^2(t)  \right\rangle \rightarrow R_0^2\;(t/\Delta t)^{2-\gamma} .

For example, directional correlations decaying with lagtime as an inverse square root will produce a MSD growing with a fractional powerlaw exponent of 1.5.

Stochastic Process Generator with arbitrary PDF and ACF

2009-07-15T17:02:00.008+02:00

Such an algorithm would be very useful for many purposes. Some ideas on this problem have been published by Primak.

Here some inefficient, but general concept:

First, a large number of independent sample values x_n (n=1...N) are generated, in accordance with the prescribed probability distribution function ( PDF, P(x) ). If the x_n would be stringed together in their original, independently generated sequence, the autocorrelation function ( ACF, Cxx(t) ) would be that of white noise. However, if the N sample values are re-shuffled in the right way (i.e. if a permutation is applied to the order of values), the ACF will change, while the PDF remains constant. It is therefore possible to use an evolutionary optimization algorithm for the re-shuffling procedure, which tries to mutate the order of the sample values until the ACF matches best the goal Cxx(t). The elementary mutations could, for example, consist of the exchange of compact intervalls of sample values.

A further Google-search yields the following, related papers:

Generation of pseudo random processes with given marginal distribution and autocorrelation function ( A. S. Rodionov et al. (2009))
Windowing to simultaneously achieve arbitrary autocorrelation characteristics and probability densities in noise generators ( R. Taori et al. )

New paper accepted for publication in PRE

2009-07-15T10:02:00.001+02:00

Title: "Noise and critical phenomena in biochemical signaling cycles at small molecule numbers"

Authors: C. Metzner , M. Sajitz-Hermstein , M. Schmidberger , B. Fabry

Abstract: Biochemical reaction networks in living cells usually involve reversible covalent modification of signaling molecules, such as protein phosphorylation. Under conditions of small molecule numbers, as is frequently the case in living cells, mass action theory fails to describe the dynamics of such systems. Instead, the biochemical reactions must be treated as stochastic processes that intrinsically generate concentration fluctuations of the chemicals. We investigate the stochastic reaction kinetics of covalent modification cycles (CMCs) by analytical modeling and numerically exact Monte-Carlo simulation of the temporally fluctuating concentration. Depending on the parameter regime, we find for the probability density of the concentration qualitatively distinct classes of distribution functions, including power law distributions with a fractional and tunable exponent. These findings challenge the traditional view of biochemical control networks as deterministic computational systems and suggest that CMCs in cells can function as versatile and tunable noise generators.

Here is the PDF of the final version on arxiv.

Talk: Noise in biochemical signaling cycles - Is life more like classical music or like free jazz ?

2009-07-13T14:37:00.002+02:00

Here is the PDF.

John Holland: Hidden Order

2009-07-04T16:07:00.011+02:00

I started to read John H. Holland's book "Hidden Order - How Adaption Builds Complexity". A few major ideas of this work will be summarized in the following post.

Cities, immune systems, central nervous systems and ecosystems as examples of Complex Adaptive Systems (CAS). CAS consist of interacting agents. Each agent is described by stimulus-response rules. Clusters of rules can generate any behaviour that can be computationally described. The range of available stimuli and responses determines the set of possible rules. Agents adapt to their environment by changing their stimulus-response rules based on experience: Learning. A large part of the environment consists of other agents.

There are 7 basic properties or mechanisms common to all CAS:

Aggregation

In the first sense: In order to model a system, we aggregate similar things into categories and then treat them as equivalent. In the second sense: Agents aggregate to meta-agents with new emergent properties, giving rise to hierarchical structures typical of CAS.

Tagging

Often several agents in a system appear indistinguishable from the outside, i.e. without closer interaction their differences are hidden. The CAS can break this symmetry of the semi-identical objects by tagging them (e.g., writing symbols on identical looking white balls). Also, different looking objects can be tagged as identical for other purposes. Tagging facilitates selective interactions and thus provides a basis for filtering, specialization, cooperation and hierarchical structure formation.

Nonlinearity

Using standard definition of linear function f(x,y,z..) = ax + by + cz +... Example for nonlinear interaction: Bimolecular reaction rates, such as f(x,y) = d/dt z(t) = a x(t)y(t). Nonlinearities typically prevent mathematical aggregation (model simplification by coarse-graining). This section of Holland's book is slightly confusing, in my opinion. In order to demonstrate the impossibility of coarse-graining a nonlinear system in the same way as a corresponding linear system, I would prefer an example such as this.

Flows

Typically, a CAS has a network structure of nodes (the agents) and links (the selective connections between agents). The dynamics of the CAS consists, beside the occasional addition or removal of nodes and links, mainly in a flow of some "material" (bio mass, chemicals, signals, goods, money, ..) through this network. The agents process the incoming flux of material into new forms and distribute it to their output nodes. Holland mentiones two general effects connected to flows. First, the "multiplier effect": If one injects additional material into a specific node, this node will typically retain a part of the extra amount and pass on a fraction r to its output nodes. In a linear chain of identical agents, the total change in the system created by the injection would be 1+r+r^2+r^3+... = 1/(1-r) > 1. So the initial change is multiplied. Second, the "recycling effect": If extra material is injected into a node that is part of a (sub-) network with loops, then the circling of the material can further increase the total effect dramatically. Also, due to circling the material is remaining longer within the system (Example: Food webs in rain forests).

Diversity

Diversity, i.e. the emergence of ever different types of agents and interactions, comes about because the agents actively try to find and occupy new niches by adaption (mutation and selection). After removing an agent, such niches are automatically refilled by similar agents (persistent patterns). With each new type of agent, further opportunities for interactions are created. In systems with recycling the number of opportunities for participation are particularly large.

Internal Models

Complex agents learn to behave in a given environment in a way that increases the probability of future reward. This requires to build an internal model of (relevant aspects of) the world, allowing to predict the results of actions. Models can be tacit (hardwired into the agent's stimulus-response machinery) or overt (explicit internal comparison of alternatives). They can be improved by adaption (variation and selection).

Building Blocks

CAS in complex environments have to respond to ever changing situations. They deal with this variety by decomposing novel situations into simpler, re-occuring and thus reusable building blocks. With a modest repertoire of nA variants of type A building blocks and nB variants of type B building blocks one can already compose nA*nB different composed objects, so this method is extremely efficient.

Emergence of long-time correlations in reaction networks - A concept

2009-05-27T11:30:00.004+02:00

I am posting here some preliminary ideas on the design of a biochemical reaction network that produces concentration pulse trains with a powerlaw distribution of pulse widths. I don't aim to find the most elegant solution, a proof of principle would be ok.

The idea is to utilize the intrinsic randomness of single-molecule reactions and to amplify it into macroscopic population fluctuations using nonlinear effects. The concept goes as follows:

There is a huge reservoir of some substance X0. It can be converted into another molecular species X by the help of an (activated) growth factor G*. Another enzyme, called the decay factor D, reconverts X back into X0. The growth reaction is assumed to be autocatalytic:

X0 + X + G* --> 2X + G*

For the decay reaction we assume that the enzyme D is operating within the saturation regime, i.e. at constant conversion rate:

X + D --> X0 + D.

At the beginning of each pulse, the autocatalytic growth reaction produces an exponential "concentration explosion", as long as the activated growth factor G* is present. If we assume, in the extreme case, that we have only a single G*-molecule around which looses its activation status spontaneously (Poisson process)

G* --> G,

we have the exponential growth of X terminated after an exponentially distributed time interval DeltaT_grow. This creates a power-law-distributed maximum concentration of X. After this event, the decay reaction reduces the concentration X at constant rate. The time period DeltaT_decay until the concentration is back to "ground level" is, therefore, also power-law-distributed. This is the main mechanism to create long-time correlations.

(The decay reaction is already proceeding during the exponential growth phase. But this should be no problem as soon as the growth rate greatly exceeds the decay rate)

Since we want to create a more or less binary concentration pulse, we couple the heavily fluctuating X-concentration to another chemical species E:

X + E <---> XE

Assume that the total number #E_tot of E-molecules is small. Then, as long as the number of available X-molecules is larger than a certain threshold, basically all E will be bound in XE-complexes, i.e. the number of XEs saturates at the limit #E_tot. Vice versa, if the X are below threshold, the number of XEs falls to zero. So the XE-signal has the desired binary pulse character. The "on-time" of XE(t) should be power-law-distributed.

There remains (at least) one problem with this scheme: How is the next pulse initiated ? To do this, we only need to re-activate the growth factor. However, it is important that this does not happen before the last pulse is over, i.e. before X is decayed back to ground level. We therefore couple the activation to the binary switch molecule E/XE:

G + E --> G* + E

E is not available during the period when X is large, because all E is bound into XE-complexes.

Cell invasion (6) : First test of chemotactic model

2009-05-25T15:27:00.025+02:00

I have implemented the coupled drift-diffusion equations using the simple numerical method tested recently.

For the first tests (program: invade3), the general parameter settings were as follows:

len=100.0;
dz=1.0;
timeSim=1000.0;
TMX=5;

As the initial conditions, the density profile of the guiding substance is set zero g(z,t=0) = 0. The cell density profile c(z,t=0) is set to a half-Gaussian of variance=1, with the maximum at the left boundary.

Run A: Nonmoving, insensitive cells. Effect of guide production, diffusion and decay

The cells remain in their initial distribution (because of D_c=0), sharply localized at the left boundary.

The guiding substance is constantly produced by the stationary cells (almost a point source). It also diffuses and decays. This leads quickly to a stationary, exponential concentration profile:

Run B: Additional slow cell diffusion, no sensitivity

Keeping all other parameters constant, the cell diffusion constant is next set to D_c=0.1. The cells, being not sensitive to the guide, now show the expected Gaussian diffusion profiles:

The initially exponential profile of the guiding substance is now gradually transforming into a Gaussian profile. This is due to the finite memory of the guide distribution and its constant re-production by the cells, which themselves assume a Gaussian distribution:

Correction: In the above figure, the green and blue curve corresponds to t=200 and t=400, respectively.

Run C: Adding sensitivity to the cells

Keeping all other parameters as in run B, we now set the sensitivity parameter to sens=3. Interestingly, the cells, starting from the initial Gaussian, now develop an approximately exponential profile:

The distribution of the guiding substance develops a non-Gaussian tail as well:

Correction: In the above figure, the lowest two curves correspond to t=200 and t=400, respectively.

The preliminary interpretation is as follows: Before the slow cell diffusion is setting in, the faster dynamics of the guiding substance is developing a concentration gradient. This gradient causes an effective back drift of the cells.

Conclusion: The chemotaxis model can produce approximately exponential invasion profiles, as observed in the experiments.

Idea for analytical case:

In the limit of D_g being much larger than D_c, the stationary g(z) should become close to a perfect exponential with a large spatial decay constant. If this decay constant exceeds the width of the stationary cell profile, it can (for those small z) be Taylor-approximated by a linear profile. The gradient of a linear profile is constant, corresponding to backward drift with constant velocity. This can be shown analytically to result in an exponential cell density profile.

It also seems mathematically reasonable that an exponential ansatz for both g(z) and c(z) can self-consistently solve the stationary coupled drift-diffusion equations (containing only derivatives and products of g(z) and c(z)) in the limit D_g>>D_c.

Note that this limit is also biologically reasonable, because small signalling molecules can diffuse much more easily through a dense matrix than the huge living cells.

Cell invasion (5) : Cooperation via Chemotaxis

2009-05-25T09:49:00.011+02:00

Back

It is known that cells can communicate with each other via signalling molecules that are diffusing in the extracellular matrix (e.g. cytokines). Would such chemotaxis be compatible with the apparent backward drift component observed in the motion of invading cells ?

To include chemotaxis into the cell invasion model, we assume that the effective drift velocity is proportional to the local gradient of the "cytokine" concentration. This leads to the following set of coupled equations:

Drift-diffusion of cells:

\frac{\partial}{\partial t}c(z,t) = - \frac{\partial}{\partial z}\left[v(z,t)\;c(z,t)\right] + D_c \frac{\partial^2}{\partial z^2}c(z,t),

Gradient controlled drift velocity:

v(z,t) = s\; \frac{\partial}{\partial z} g(z,t)

Production, diffusion and decay of the guiding substance:

\frac{\partial}{\partial t}g(z,t) = p\!\cdot\!c(z,t) + D_g \frac{\partial^2}{\partial z^2}g(z,t) - d\!\cdot\!g(z,t)

Here, the two coupled fields are the cell concentration profile c(z,t) and the concentration profile of the guiding substance (the "cytokine") g(z,t). The diffusion constants for cells and guiding substance molecules are denoted by D_c and D_g, respectively. In addition, we introduce a sensitivity parameter s, a production parameter p and a decay parameter d. It is assumed that the guiding substance is locally produced by each cell at a constant rate and that it decays within a characteristic time 1/d.

Cell invasion (4) : Numerical solution of drift-diffusion equation

2009-05-22T15:06:00.017+02:00

Back

The cell invasion problem requires an accurate numerical solution of coupled partial differential equations. The fields (concentration profiles) f(z,t) will depend on one spatial coordinate z and the time t. The spatial coordinates will be restricted to the intervall [0,len], with well-defined boundary conditions at the left and right limits: The current must dissapear outside the intervall.

I will try to solve this problem by discretizing the z-coordinate into a regular grid,


z_k = k\;\Delta z,

and then treating each z_k as a separate dynamical variable. The problem can thereby be mapped onto a set of coupled ordinary differential equations, which are solved with the Runga-Kutta method.

As a simple test and practice case, consider the drift-diffusion equation


\frac{\partial}{\partial t}c(z,t) = - v \frac{\partial}{\partial z}c(z,t) + D \frac{\partial^2}{\partial z^2}c(z,t),

with position- and time-independent diffusion constant D and drift velocity v.

For the following tests, we set the initial distribution c(z,t) equal to a narrow, normalized Gaussian (variance=10), placed at the center of the intervall. The spatial resolution is set to dz=1.

Test 1: Free diffusion in a box

First we set the drift-velocity equal zero. The initial peak should broaden with time. Asymptotically, a flat distribution should develop, while the integral (total particle number) is preserved. The simulation shows the expected behaviour:

Test 2: Diffusion with drift (boundaries far away)

The case of combined diffusion and drift has been compared directly with the analytical solution (solid lines). One expects the center z_c(t) of the Gaussian to move according to


z_c(t) = v\;t

and the variance to grow linearly


\sigma^2(t)=\sigma^2(0)+ 2Dt

The simulation (dashed lines) fits the analytical curves nicely, even for the modest spatial resolution of dz=1:

Test 3: Diffusion with drift against a wall

Finally, we choose a negative drift velocity and wait until the distribution is pushed againts the left wall. We now expect the stationary solution


c(z)=\frac{v}{D}\;e^{-(v/D)z}.

Again, the simulation (symbols) agrees very well with the analytical solution (solid lines):

As a conclusion, the simple numerical method seems to work as expected. I have also checked that the results do not change significantly when the resolution dz is reduced.

Next

Cell invasion (3) : Effect of cell proliferation

2009-05-07T15:27:00.011+02:00

Our biologists tell me that there might be some cell division be going on during the invasion process. I included a corresponding proliferation term to the master equation.


P_z^{(t+1)}=\;\ldots\; +\alpha P_z^{(t)}.

In the figures below, dashed lines correspond the case with zero growth rate, solid lines to a growth rate of alpha=0.01.

----------------------------------------------------

Symmetric case

----------------------------------------------------

Drift to the right (chemo-attractant)

----------------------------------------------------

Drift to the left (chemo-repellant)

----------------------------------------------------

Generally, the effects on the z-profiles are more of a quantitative nature. Drastic differences result for the surface cell density versus time.

Cell invasion (2) : First experiments

2009-04-30T15:30:00.012+02:00

Experiments on the trans-endothelial migration and collagen invasion of tumor cells have been published in a paper by C. Mierke et al. in Biophys J. 94, 2832 (2008).

With their assay,

the authors meassured the invasion depth of the cells into a 3D collagen gel fiber network. They compared two cases. In case 1 (dark gray bars), the cells first had to transmigrate through a HUVEC endothelial barrier before they could invade into the collagen. In case 2 (light gray bars), the barrier was missing. The results are shown below:

A semilog-plot of the data suggests an exponential invasion profile:

This behaviour would correspond to the "chemo-repellant" case of the model posted as Cell invasion (1).

Cell invasion (1) : Biased diffusion model

2009-04-29T16:01:00.033+02:00

Consider living cells that are initially located at the surface (x-,y-plane) of a half-space of dense medium (positive z-axis of coord. system pointing into the medium, z=0 at the surface). As shown in the figure below, the cells can invade into the volume of the medium, migrate there and possibly return to the surface.

As the simplest possible model, we assume that the cells perform a (biased) random walk with statistically independent spatial components (dx,dy,dz) of the displacements.
We are only interested in the temporal evolution of the probability density function P(z,t) of the cell's z-position. Since the 3 coordinates are independent random processes, we consider the whole problem as a 1D-diffusion problem with special boundary conditions and with possible trends (drift terms).
The trends could, for example, describe the effect of a chemo-attractant in the medium, creating a bias of diffusion into the positive z-direction (later called the "right" direction). A chemo-repellant would cause a drift to the "left".

The corresponding Fokker-Planck equation for the temporal evolution of P(z,t) could be solved analytically. Instead we discretize the problem spatially (z-segments in the figure) and temporally and integrate the resulting master equation numerically. The surface is described by the special segment Z=0.

For given resolutions (dz, dt), we define 4 dimensionless parameters:

muR = fraction of particles in segment Z that move to the segment Z+1 to the right within one time step.

muL=fraction moving to the left.

muV=fraction of particles on the surface (Z=0) that move into the volume.

muS=fraction of particles in volume segment (Z=1) that move to the surface.

Let


P_z^{(t)}

denote the normalized probability density (PDF) of particles in segment z at time step t.

The initial condition is


P_z^{(0)}=\delta_{z0}.

The master equation reads for the surface segment


P_0^{(t+1)}=(1-\mu_V) P_0^{(t)} + \mu_S P_1^{(t)},

for the volume segment Z=1


P_1^{(t+1)}=(1-\mu_S-\mu_R) P_1^{(t)} + \mu_V P_0^{(t)} + \mu_L P_2^{(t)},

and for all other segments Z>1


P_z^{(t+1)}=(1-\mu_L-\mu_R) P_z^{(t)} + \mu_R P_{z-1}^{(t)} + \mu_L P_{z+1}^{(t)}.

First, we shall treat the passage of the cell through the s pecial boundary between the last segment Z=1 and the surface segment Z=0 in the same way as for all other segments:


\mu_V=\mu_R \; , \; \mu_S=\mu_L.

---------------------------------------------------------

Simulation results for symmetric conditions:


\mu_R=\mu_L=0.1

Here, the density profiles are approximately Half-Gaussians, peaked at the surface, with temporally increasing variance. The surface density decays at long times according to a t^(-1/2) powerlaw.

---------------------------------------------------------

Simulation results with "chemo-attractant":


\mu_L=0.1 \;,\; \mu_R=0.2

Here, the peak of the Gaussian profiles is drifting to the right. The surface density decays exponentially.

---------------------------------------------------------

Simulation results with "chemo-repellant":


\mu_L=0.2 \;,\; \mu_R=0.1

Here, the profiles quickly approach a stationary, exponential distribution. The surface density falls to a stationary, finite value.

---------------------------------------------------------

We will next compare those simulations with measured cell invasion data.

Confusion about nonlinearities in biochemical networks

2009-04-23T10:16:00.030+02:00

In a recent discussion I was temporally confused about the question of whether the covalent modification network

is linear or not.

(I admit that this already has made me think about the possibility of nonlinear behavior emerging as an artifact of a coarse-grained description of a fundamentally linear dynamic system.)

In this post I want to clarify this point.

The following figure shows an arbitrary example of a biochemical reaction network:

The letters stand for chemical species and at the same time for the dynamical variables describing the temporally changing concentrations of these species. Here, X is an "input chemical", i.e. its concentration is completely under control of the experimenter. The concentration of the "output chemical" Y is measured by the experimenter, without disturbing its value. Inside the reaction network there can be an arbitrary number of other species z_n, dynamically coupled by chemical reactions. Note that we describe the system, from the outside, by only one input and one output variable, although the internal dynamical state has more dimensions.

We consider y(t) as a general functional of the input x(t):


y(t) = F\left\{x(t)\right\}.

If the reaction network would be a linear system, the functional could be written as a convolution with a Greensfunction:


y(t) = \int_{-\infty}^t G(t-t^{\prime}) x(t^{\prime}) dt^{\prime}.

In a mass action approximation, the temporal change of the various chemicals is described by coupled first order differential equations:


\frac{d}{dt}y = f_y(x,y,\vec{z})

and


\frac{d}{dt}\vec{z} = f_z(x,y,\vec{z}).

Here, the vector z stands for all the interior chemicals of the network.

The only case where we can actually define a Greensfunction is the trivial case where the right hand sides of the above equations are simple linear combinations of the variables:


f_y(x,y,\vec{z})=a\;x+b\;y+\vec{c}\;\vec{z}


f_z(x,y,\vec{z})=d\;x+e\;y+\vec{g}\;\vec{z}.

Here, the quantities a-g stand for fixed coefficients. More complicated terms such as products of variables are not permitted.

What does this mean for the structure of biochemical reaction networks ? We must not have any bi- (or higher-) molecular reactions involved. For example, in the above example network, the equation for the temporal change of Z2 reads


\frac{d}{dt}\vec{z_2} = k_1\; x\; z_1 - k_2\; z_2.

The " x * z1 " term already makes the system nonlinear.

The covalent modification network also contains bi-molecular reactions (where the enzyme-substrate complexes are formed). It is therefore nonlinear.

In conclusion: There are no multi-molecular reactions in a linear network.

Adaptive fiber networks

2009-04-22T17:45:00.001+02:00

Consider a model cell in the "spider" geometry: A central cell body with N radial stress fibers, each fiber starting at the cell center and terminating at the remote end in a focal adhesion contact to the substrate. (Compare figure).

Let us assume that the fibers can be modelled as linear springs of stiffness k, however with a time-dependent rest length L0(t). The force is at any time given by F=k*(L-L0), if the actual fiber length L exceeds the rest length L0. The force is zero otherwise (stress fibers don't withstand much compression).

Now assume the fibers are adaptive. They have a preferred (or "goal") prestress F_goal, and they change their rest length in order to come closer to this goal. For simplicity, assume a constant adaption velocity in both directions:

dL0/dt = +v if F>F_goal (lengthening)
dL0/dt = -v if F<F_goal (shrinking)

Birth of fibers: The fibers are born with fixed initial length L(t=0) and with random directions. There are available "slots" for only NMX fibers in total. The birth rate of the fibers is proportional to (N-NMX) and therefore self-limiting. At time of birth, a fiber is relaxed, i.e. L0=L(t=0).

Death of fibers: Assume that fibers die if their length falls below L_min (under-length) or exceeds L_max (over-stretch). Dead fibers immediately dissappear.

What would be the expected behaviour of that model ?

Remove all fibers from the model cell and place it onto a flat substrate with friction. After some time, a first new-born fiber will appear (say, at the right side of the cell) and adhere to the substrate with its remote end. Being initially relaxed, it will start to shrink. This generates a traction force on the cell, which is thus pulled to the right side against the friction force.

Let's say that before the focal end point is reached, another fiber is born at the left side of the cell. It also starts to shrink, making it harder for the right fiber to increase tension and to reach the goal prestress. So, both fibers work against each other. The cell will now be the node of a spring network with time-dependent prestress and perform some "complicated" motion.

If not disturbed by further birth processes, the cell might reach an euqilibrium position, with both fibers in their goal prestress state.

What if the second fiber also appears at the right cell side ? Then both fibers shrink until they die from under-length. Fibers need some resistence to develop their goal prestress. This resistance can be provided by other fibers or by some external forces (if the cell is sitting in a mechanical potential well).

Consider next an initial situation with N fibers, all in their ideal state (F=F_goal).

If we quickly displace the cell center from its equilibrium position, we shall experience spring-like restoring forces from those fibers that are stretched, zero forces from those that are compressed. All in all an elastic response (With some viscous contributions).

If the displaced cell is quickly released, and if the displacement was not so large to cause the death on any fiber, the cell will return to its equilibrium postition because there was not enough time for adaption.

But if we hold the cell for a long time in the displaced position, the fibers will adapt their rest lengths: The stretched ones will lengthen, the compressed ones will shorten. As a result, after releasing the cell it will return to a new place in between the former equilibrium position and the externally imposed position (plastic deformation).

Project week: Anja Michl

2009-04-10T10:47:00.003+02:00

A new one-week project has started with Anja Michl. We plan to study Self-propelled agents in complex potential landscapes. The running project will be reported (in German) in a dedicated research blog.

Scale-free distribution of molecule numbers in signaling cycles

2009-04-06T16:40:00.002+02:00

I have written a new paper manuscript on the following topic:

Biochemical reaction networks in living cells usually involve reversible covalent modification of signaling molecules, such as protein phosphorylation. Under the frequent conditions of small molecule numbers, mass action theory becomes insufficient to describe the dynamics of such systems. Instead, the biochemical reactions must be treated as stochastic processes, producing intrinsic concentration fluctuations of the chemicals.

We investigate the stochastic reaction kinetics of covalent modification cycles (CMCs) by analytical modelling and numerically exact Monte-Carlo simulation of the temporaly fluctuating concentration x(t). The statistical behaviour of this simple network module turns out to be so rich that CMCs can be viewed as versatile and tunable noise generators. Depending on the parameter regime, we find for the probability density P(x) several qualitatively different classes of distribution functions, including powerlaw distributions with a fractional and tunable exponent. These findings challenge the traditional view of biochemical control networks as deterministic computational systems.

Bacterial chemotaxis as a guiding line

2009-03-02T14:36:00.014+01:00

We still lack a coherent, integrative understanding of how eukaryotes sense the properties of their environment, how they process this information, and how exactly they respond in re-organizing their cytoskeleton and thereby move in a goal-oriented way.

In contrast to that, our understanding of bacterial chemotaxis is much further advanced. We know the main players of the chemotaxis signal network, from the aspartate receptor-complex downstream to the flagella motor. One can actually compute, in a quantitative model, how a change in the extra-cellular concentration of attractant molecules eventually affects the relative fractions of clockwise and counter-clockwise rotations of the motor and thus leads to a change in the bacteria's tumble frequency. We understand how this modulated random walk brings the organism closer to the attractant.

People have also successfully simulated the effect of knocking out certain proteins in the signal pathway. One even understands how robust these biochemical signal processors are with respect to random parameter variations.

Why is it so far not possible to achieve a similar understanding of eukaryotic cell migration ? Why can't we take bacterial chemotaxis as a guiding line for developing a model ?

There are so many obvious similarities.

For example, the whole strategy for approaching the goal (attractant) is based on an internally generated random process (stochastic tumbling events in between segments of swimming along a straight line), which is only modulated by the signals of the receptor. (Note that this is in line with the recent insights how stochasticity is exploited in biochemical systems.) This search strategy seems related to the exploratory behaviour that I suggested for migrating cells in an environment with sparse adhesion opportunities and hindrances.

It has even been found that there are long-time correlations in the statistics of the tumbling intervalls. The authors have mentioned the possibility that the resulting Levy-walks might reflect an optimized search strategy.

Subsequently, other authors have tried to indentify a biochemical mechanism that might produce the necessary power-law correlations in the motor control signal.

Converting two short-time-correlated signals with powerlaw PDF into one long-time-correlated signal

2009-02-27T14:13:00.040+01:00

Assume a stationary, non-negative signal u(t)>=0, exponentially autocorrelated with characteristic time constant \tau_c


C_{uu}(\tau)=\left< \left(u(0)-\overline{u}\left) \left(u(\tau)-\overline{u}\left) \right> = \sigma_u^2 \;e^{-\tau/\tau_c},

having a powerlaw-shaped probability density function (PDF) with exponent a in the range [2,3], up to some upper cutoff value u_max:


P(u) \propto \theta(u_{max}\!-\!u)\;u^{-a}.

The cutoff is necessary for a finite variance of the PDF.

Assume a second random signal d(t), independent or u(t), with the same statistical properties.

From u(t) and d(t) we generate a new signal y(t) in the following way:

y(t)=\int_0^t \left[u(t^{\prime})-d(t^{\prime})\right] dt^{\prime} = \int_0^t \left[x(t^{\prime})] dt^{\prime}.

Note that x(t) is a random signal with positive and negative values, symetrically distributed around zero, with powerlaw tails. What is the autocorrelation function of y(t) ?

To build up some intuition, we can approximately replace x(t) by a simpler function z(t) which is piecewise constant over periods of lengths \tau_c:


z(t)=\sum_n \theta(t-n\tau_c) \;\theta((n\!+\!1)\tau_c-t)\;x_n,

where the heights of the steps are defined by


x_n = x(n\tau_c).

When this z(t) is integrated (instead of x(t)), the output signal y(t) will increase linearly within each step. At the end of each step, i.e. at discrete times


t_n=n \tau_c,

z(t) has performed a step of width


\Delta y_n=x_n\tau_c,

which is, of course, powerlaw distributed. In other words, embedded into the random process y(t) there is a random walk with constant waiting times and powerlaw-distributed step lengths - a Levy flight:


P(\Delta y) \propto \Delta y^{-a}.

It is known that such a Levy flight has a mean squared deviation (MSD) of


\left< \left( y(t)-y(0) \right)^2\right> \propto t^{2/(a-1)}.

So, if a=3, then the MSD of y(t) behaves diffusive (Brownian walk). If a=2, then the y-MSD behaves ballistic. For values in between this range, one obtaines a signal y(t) with fractional, superdiffusive behaviour. This, in turn, means that y(t) is long-time-correlated.

The above was only a hand-waving argument, which needs to be tested.

But if it turnes out to be true, it provides a simple biochemical mechanism for producing a long-time correlated fiber growth process: As an input, one only needs a reaction network that produces fluctuating concentrations u(t)/d(t) of proteins U/D, with statistical properties as the u(t) above (We actually know examples of such networks !).

The protein U acts as an enzyme for catalyzing the assembly of the stress fibers S from a huge reservoire of raw material R. Protein D acts a a disassembler:

R + U ---> (R-1) + S + U
(R-1) + S + D ---> R + D

So the momentary growth rate of the fiber will be proportional to x(t). Then, the signal y(t) above can be interpreted as the fiber size. We have already proposed a similar mechanism in Long Beach, California, 2008 as our "Model 3".

Maybe it is possible to replace the fluctuating disassembly process by one occuring at constant rate, such as described in post IMCM[3].

Powerlaw Foraging

2009-02-27T10:20:00.008+01:00

Many foraging animals try to find their prey by random walks, i.e. chains of straight segments (steps) separated by random turns. As measurements show, these walks are not Brownian, but correspond to Levy flights with powerlaw-distributed lenghts L of the straight segments:


P(L)\propto L^{-m}.

Here, the exponents m range from 1 to 3 and the resulting stationary PDF of visited sites P(x,y) is not Gaussian.

In a 1999 Nature paper, G.M. Viswanathan shows that such powerlaw foraging (with exponent m=2) can be an optimum search strategy if target sites are sparse and can be visited repeatedly.

For an intuitive understanding, he notes:

In a Levi flight, irrespective of the exponent, the probability of returning to a previously visited site is smaller than for a Brownian walk with Gaussian PDF.

Also, N Levy walkers visit much more new sites than N Brownian walkers.

IMCM [4]: Why long-time correlated fiber growth rates ?

2009-02-25T17:22:00.007+01:00

BACK

Here my very hand-waving argument:

When R_ass(t) is long-time correlated, the fiber strength s(t), which is essentially an integral over R_ass(t), has also a power-law autocorrelation function. Since the momentary cell position depends linearly on the strengths s_j(t) of the individual fibers, the cell as a whole will perform a persistent, power-law random walk (i.e. the mean square displacement increases with lagtime as a fractional powerlaw).

So what ?

First of all, a fractional persistent random walk has actually been observed for migrating cells.

Second, it has been shown in foraging theory that such non-Brownian walks are optimal search strategies when the targets are sparsely distributed.

Third, there are simple biochemical reaction networks that produce long-time correlated growth rates R_ass(t).

However, I don't want to insist too much on powerlaws as the optimum strategy. Computer experiments could be used to systematically compare the efficiency of the resulting cell migration in given environments when different statistics for R_ass(t) are used.

IMCM [3]: The life of stress fibers

2009-02-25T16:19:00.008+01:00

BACK

Let us define the birth of a stress fiber as the moment when a protrusion with preliminary adhesions has "passed the traction test" and starts to grow.

In order to have a clear working model to start with, we can assume that the fiber is subject to disassembly processes at a fixed rate R_dis and to re-assembly processes at a variable rate R_ass(t), the latter depending on various factors (to be discussed later). The "size" or "strength" s(t) of the fiber shall change according to

d/dt s(t) = R_ass(t) - R_dis.

The fiber is alive as long as s(t) exceeds some minimum size s_min. When R_ass(t) is smaller than R_dis for a too long period, s(t) falls below s_min and the fiber dies. After this, its remaining material is recycled by the cell and the fiber "slot" becomes available for a new protrusion. If the assembly rate R_ass(t) has some randomly fluctuating component, the resulting life time of each fiber becomes a random variable, too.

If we now remember our main working hypothesis that the cell has the chief goal (is evolutionary optimized) to migrate efficiently in the complex environment, we come to an interesting question:

How should the assembly rate of each fiber j be controlled, or regulated, in order to optimize migration ?

Is it neccessary for the cell to couple the dynamics of different fibers, analogous to the coupled control of legs in a spider ?

The most simple and biologically un-expensive way of "control" would be that the assembly rates of the individual fibers j are mutually independent, randomely fluctuating functions R_ass_j(t). In this case, the question is: Which statistical properties (distribution function, autocorrelation function, etc.) are optimum for efficient migration in environments of given statistical properties.

For reasons that I shall discuss later, I have the hypothesis that R_ass_j(t) with long-time-correlations (i.e. the autocorrelation function decays with delay time as a power law with fractional exponent e) are the optimum choice for the cell.

NEXT

IMCM [2]: Exploration of environment, selection of adhesions, chemotaxis

2009-02-25T10:23:00.007+01:00

BACK

IMCM starts at top level with a clear goal for the cell: To find a passable migration path within the complex environment. This environment contains opportunities (HSs) as well as hindrances (EOs), and it requires some minimal intelligence to find a solution to this problem.

In such situations, nature often reverts to the evolutionary trial-and-error method: Create a random variety of proposals, evaluate them, then keep the good and reject the bad ones.

Within the cytoskeleton, to provide a known example, the spatial arrangement of microtubuli (MT) is not static, but subject to ongoing dynamic reorganizations. Even if the overall MT arrangement appears constant over extended periods, the individual microtubuli are building and disintegrating on a shorter time scale. They are growing radially outward from some nucleation center, in random directions. But only those individuals that find a stable anchoring point at the other side are stabilized. The others disassemble more quickly. So the seemingly constant arrangement is just a stationary dynamic state. The same exploratory behavior is found again and again in biology and other self-organizing systems.

It therefore seems reasonable to assume a similar mechanism for the cell migration problem: Protrusions into random directions are "seeking" for hard substrate patches.

During the exploratory phase of a PT, preliminary adhesions are formed with the substrate. Next, the PT is applying increasing contractile forces to the adhesions, in order to evaluate the stiffness of the patches.

In the case of soft substrate, no significant tension will build up in the PT, and the unbalanced contractile force will simply pull back the PT to the cell body. But on hard substrate, there will develop tension, and this tension can be used as a signal to reinforce the traction strength of the PT. Only such "successful" PTs enter into the mature phase and become stress fibers with a possibly much longer life time.

Clearly, this scheme requires a kind of molecular tension-sensor, located in the PT-fiber or in the adhesion. Under tension, the sensor changes conformation and activates the biochemical reactions leading to reinforcement of the fibers.

The exploratory mechanism described so far allows the cell to find (even sparsely distributed) hard patches in the close environment and to build stress fibers, firmly connecting the cell body with those anchoring patches. In the model, we can also allow for multiple fibers connecting to the same hard substrate patch (One often observes parallel bundles of fibers).

However, it is reasonable to assume a maximum number of PT/SF-"slots" for each cell body. So let the rate of formation of new exploratory PTs be proportional to the number of currently available "slots".

If the mature stress fibers would live forever, the cell would end up in a state similar to a tree, firmly rooted to the ground, but not mobile. The existing fibers can still fluctuate in their traction forces, leading to bounded spatial fluctuations of the cell body.

But true cell migration requires that stress fibers have a finite life time and can be replaced by new ones, pulling the cell into other directions. The long time dynamics of mature fibers will be discussed in another post.

Even if the cell can migrate efficiently, thanks to an optimized life time statistics of the mature fibers, the (statistically averaged) migration will still be radially symmetric in a homogeneous and isotropic environment.

In such environments it may be necessary for the cell to have a vague sense of the direction. Let us therefore keep in mind the possibility of some chemotactic mechanism. For the simulation, we could simply assume that the probability of formation of new exploratory PTs is slightly higher in the directions of the chemo-attractant.

Finally, I would like to point out that the above exploratory mechanism needs not to be simulated in detail. For the migration process, the only thing that counts are the surviving, stable stress fibers. Thus, it may be sufficient to ascribe an effective rate for stable stress fiber formation to any of the closeby hard patches (Small problem: They can be shadowed by obstacles), with a global rate proportional to the cell bodie's available number of slots.

NEXT

IMCM [1]: Cell in a complex environment

2009-02-24T17:16:00.011+01:00

BACK

I start with a 2D model, but the generalization to 3D is straight forward.

Let the main body of the cell be an elastic ball (viscous properties can be added later) of a certain diameter D_cell. The cell can be squeezed in order to pass through environmental constrictions, but in the relaxed state this is its standard size. Originating from the main cell body, there can grow radial protrusions (PTs). I'll come back to that point soon.

The environment is inhomogeneous. Consider a 2D regular grid of quadratic patches. Each patch can be in any of three different states: Soft substrate (SS), hard substrate (HS), or elastic obstacle (EO).

A cell can freely move parts of its body (which is larger than a single patch) over soft SS patches.

Elastic obstacles, however, produce a local potential barrier (finite range and height, radial symmetric, centered at the respective EO patch). The potential wells of the many individual, randomely distributed EOs add up to a complex potential landscape with tunable parameters. Mathematically, the force between the cell and each individual obstacle is simply a function of the distance between cell center and EO center, and all EO forces add up.

The cell can form stable adhesions only at the hard substrate patches. When a radial cell protrusion finds a HS patch, it forms an initial adhesion there. This stabilizes the protrusion and it has then a chance to develop into a stress fiber (SF).

Some final remarks: It may be advantageous to describe the repulsive EOs by potentials that drop (sharply yet smoothly) at the borders of the obstacle, but continue to decay on a low level until infinity. That way, the total potential landscape of all obstacles is a continuous field and has well-defined minima (in the absence of stress fibers), where the cell will naturally rest.

The formation of stress fibers can be viewed as the self-creation of attractive potentials by the cell, which tend to compensate for the repulsions. So, in order to move, the cell is actively reshaping its surrounding potential landscape in a time-dependent manner.

By the way: The cell can, as another known strategy, actively reduce the height of repulsive potentials by chemically dissolving the obstacles. This effect could also be simulated rather easily.

NEXT

IMCM [0]: Proposal of a new project "Integrative Model of Cell Migration"

2009-02-24T17:11:00.010+01:00

In biophysics, systems are overwhelmingly complex. There are so many possibly relevant - yet not precisely known - interactions between the components that one cannot simply attempt to compute the "properties of the given system". There is no well-defined system to start with !

In order to ask meaningful questions, biologists have long ago learned to focus on purpose and function. Once you know the useful function a system is supposed to perform, you get some idea which components and interactions might be absolutely relevant for that purpose and which are only of secondary importance.

For this reason, I believe that also in the field of cell mechanics, an integrative, function-oriented model would be very useful, no matter how simplistic. Such a functional multi-scale model should, in particular, make a logical connection between molecular reactions, cytoskeletal reorganizations at the mesoscopic level, and the final goal of a cell: to migrate and navigate in a complex environment.

My title for this new project will be "Integrative Model of Cell Migration". The posts related to the project will be denoted by IMCM [number].

IMCM will be minimalistic. Its main purpose is to provide a concrete visual model as a basis for further discussion. It will produce testable predictions and, hopefully, induce many new questions.

My ideas about the design of IMCM will be described in the following posts.

FORWARD

[6] Project Week: Michael Schmidberger

2009-02-13T19:18:00.005+01:00

I just finished a one week theory project with Michael Schmidberger on Fluctuations in Biochemical Reaction Networks. The project has been reported online in Michael-Blog. The text is in German, so far, but I plan to post the main results here in the near future.

[5] Effect of binning and flipping

2008-12-18T14:31:00.006+01:00

The following relates to [2]-[4]:

Martin has pointed out that the natural binning of the raw data and the different flipping operation of Gonzalez et al. could be partly responsible for the marked differences between their and our results. In particular, the double-peak structure might be reduced if our raw data are artificially binned and the same flipping rule (most occupied cell at right hand side) is used.

When applied to the experimental microbead data, the results are considerably closer to the single-peaked, triangular distributions of Gonzalez:

[4] Finite Size PDF of independent points in 1D

2008-12-17T12:02:00.042+01:00

The following relates to posts [2] and [3]:

--------------------------------------------------------------------------------

Simplifications:

In order to indentify the minimum conditions of a double-peaked universal PDF, the situation can be simplified in several aspects:

A) 1D space:

By restricting the trajectories to one spatial dimension x, rotations are not required any longer.

B) Sets of N independent points:

Trajectories are normally generated by subsequently adding random increments (steps) to the respective last position of a walker. Here, in order to avoid any correlations between successive points, the N positions in each point set are drawn independently from a fixed probability density.

C) Random points equally distributed in [0,1]:

Let this probability density be constant in the intervall [0,1] and zero outside.

--------------------------------------------------------------------------------

Statistical quantities of point sets:

For each point set {x_i}, the center of mass is computed,


\overline{x} = \left\langle x_i \right\rangle_i = \frac{1}{N}\sum_{i=1}^N x_i \;,

as well as the standard deviation


\sigma = \sqrt{\left\langle (x_i-\overline{x})^2\right\rangle_i} \;.

--------------------------------------------------------------------------------

Transformations on point sets:

C-Operation (Center):


x_i \rightarrow (x_i-\overline{x}) \;\;\forall i

S-Operation (Scale):


x_i \rightarrow \overline{x}+(x_i-\overline{x})/\sigma \;\;\forall i

Note that C and S commute with each other.

--------------------------------------------------------------------------------

Effects of transformations:

The effects of C and S on the resulting averaged distributions P(x) are discussed in the following:

* No C, no S:

Direct averaging, without any C or S, yields the expected box-shaped distribution in [0,1], centered around 1/2.

* Only C:

Applying only the C-operation yields PDFs centered around zero. For N=1 one obtains a delta-function, for N=2 a triangular function, etc. For very large N, the box is recovered:

* Only S:

Scaling alone produces already a double-peak structure, centered around 1/2. However, it is a finite size effect that quickly disappears for long trajectories:

* C and S:

The combined effect of C and S yields a double-peak centered around zero:

If instead of the box-shaped distribution, a Gaussian distribution is used for the random points, the results are qualitatively the same. However, the double peak is much weaker pronounced and visible only for very small trajectory lengths.

--------------------------------------------------------------------------------

Summary:

Summing up, the double-peak is not a very remarkable feature of 1D trajectories.

--------------------------------------------------------------------------------

Origin of double-peak:

Take the extreme case of N=2: After a C-operation the two points lie symmetrically left and right from x=0. An additional S-operation scales the distance between the two points to the norm value, and the average distribution P(x) consequently consists of two delta-functions.

On the other hand, when N becomes very large, the histogram of each individual trajectory will already reflect the ensemble average very closely. The CS-operations therefore do not change the shape of this distribution qualitatively and one expects for P(x) to recover the fundamental distribution (in our case: box-shape).

Then, for intermediate lengths N, one expects a gradual interpolation between the double-delta-peak and a box-distribution.

--------------------------------------------------------------------------------

[3] Universal PDF of 2D Brownian and Persistent Trajectories

2008-12-17T08:18:00.019+01:00

I am grateful to Martin Reichelsdorfer, who has programmed and applied the CRFSA-method (center, rotate, flip, scale and average), as described in post [2], to various synthetic and experimental trajectories in 2D.

In contrast to the original method of Gonzalez et al., his raw data trajectories consisted of real points in 2D. Cells were only introduced in the final averaging step. Therefore, the flipping operation had to be redefined: It was done in such a way that the last point of the trajectory was always located right of the first point.

He first tested the method for a standard Brownian Random Walk. For a given time period dt, the PDF of a Random Walk is known to be a Gaussian, centered at the origin, with a variance proportional to dt. This is also what one naively expects to see as the ensemble averaged spatial distribution. However, after CRFSA one obtaines a double-peak structure:

Second, Martin applied CRFSA to measured trajectories of micro-beads bound to living MEVO cells. On long time scales, those trajectories show superdiffusive behaviour with fractional powerlaw exponents in the mean squared displeacement (MSD), corresponding to directional persistence in the bead motion. This persistence is a qualitative difference to the Brownian Random Walk and one might expect new features for the ensemble averaged spatial distribution function (SDF). However, the results for the Persistent Walks are not much different from before.

No significant changes are observed when the ensemble of trajectories is divided into groups of similar powerlaw exponent.

These results lead to a couple of questions:

Why is the CRFSA-SDF of a Random Walk not a Gaussian ?
What exactly is the origin of the double-peak ?
What makes the trajectories of mobile phone users so qualitatively different from Brownian or peristent walks ?

[2] Universal Spatial Distribution of Mobile Agents

2008-12-16T14:15:00.023+01:00

In a recent paper "Understanding individual human mobility patterns" published by Marta C. Gonzalez et al. in Nature, the authors tracked the coordinates of mobile phone users over an extended period.

In this study, space was divided into cells (the range of a single mobile phone antenna) and the actual raw data consisted of the sequence of cells in which each user was found. In most of the cells, the same user was found several times. In particular, there was a cell of maximum occupancy for each user.

For each individual, the spatial distribution pattern (formally equivalant to a mass distribution) was artificially centered around the origin. Next, it was rotated around its center of mass, so that the main axis of the inertia tensor became alligned along the x-axis. Next, the pattern was flipped, such that the cell of maximum occupancy was located at the right side. It was then scaled to a fixed standard deviation in both x- and y-directions. The centered, rotated, flipped and scaled distributions were finally averaged over all individual trajectories (CRFSA method) . The result was a "universal" and surprisingly non-trivial distribution pattern:

It seems worthwhile to apply a similar procedure to the trajectories of microbeads attached to the living cytoskeleton. I shall name this project "Universal PDF".

[1] Early New Year's Resolutions

2008-12-16T09:39:00.005+01:00

OK,OK - I admit that I have neglected this blog recently. Actually, I was using a Wiki system (wikidot.com) instead for documenting some of my scientific work. It has a nice Latex feature preinstalled, allows uploads of arbitrary file types, and offers more flexibility in many respects.

However, I came to see that the simple chronological ordering of posts in a blog, combined with tagging, is not so bad afterall. You don't waste time in thinking about the best way to organize your wiki sites. In addition, as the use of Latex became more convenient now and one can upload files simply from an external FTP folder (see "Online Resources" at the top right of the blog), I decided to restart blogging.

So, as a kind of early New Year's resolutions, I promise to post regularly again from now on ...

Easy LaTex-Formula in Blogger

2008-11-12T20:12:00.026+01:00

So far, I have been using the Latex Equation Editor "CodeCogs" for rendering mathematical equations in this blog. It was not very convenient. A few minutes ago I discovered "YourEquations.com" and the procedure became much easier ! For example, the HTML-code

(pre lang="eq.latex">
a^2+b^2=c^2
(/pre>

produces the formula:


a^2+b^2=c^2

Of course the "(" must be replaced by "<" in the HTML-code. Unfortunately, there is a perceptible time lag before the image of the formula appears.

Complex Systems Primer

2008-10-11T17:19:00.006+02:00

Claudius Gros from Frankfurt University has put online a preliminary version of a (still unpublished) text book on complex and adaptive systems theory. File has more than 4 MB and takes some time to download.

Stochasticity and Cell Fate

2008-10-02T15:31:00.010+02:00

A 2008 Science-Paper from Richard Losick gives some examples why stochasticity in biochemical (genetic) control networks can be usefull:

Persister state: A small fraction of a bacteria population goes into a non-growing state. So in the case of antibiotics that attack only growing cells, at least this fraction remains alive.

Swimming or chaining: Some bacteria form spatially fixed chains, feeding on local food. Another statistical group is motile and explores new food sources.

Red and green cones in human retina: A fraction of neurons become green sensitive, another fraction (non 50%, determined by distance between control region and gen) red sensitve. One fate inhibits the other. The two forms are randomely distributed over the retina.

Getting Things Done

2008-10-02T15:11:00.007+02:00

A very stimulating paper about the optimum organization of workflow:

Francis Heylighen and Clément Vidal: "Getting Things Done: The Science behind Stress-Free Productivity"

The goal of GTD is not to so much to optimally achieve given goals with given priorities (because goals and priorities are in constant flux in an information society), but rather to maximize the number of usefull tasks performed.

Cognition is

situated (triggered by environmental stimuli)
embodied (sensor- and motor-activity involved)
distributed (over the brain and environment)

The environment serves as an external memory. It triggeres actions and is a source of

affordances (opportunities)
disturbances
feedback (compare present project state with goal)

Project plans in GTD:

overall goals sufficient
not too far ahead
focus on actions

Paper-Project 1: Finite Resources

2008-07-18T15:12:00.009+02:00

Paper should contain chapter with statistics of production/consumption processes, producing Gaussian, exp-correlated concentration fluctuations
Applicable to economy (workers/construction sites), computer multitasking (CPUs/jobs), biochemistry (enzymes/assembly sites)
Check how different P(s_max) affect the P(T_compl)
Graphs: progress of indiv. jobs s_i(t), # completed jobs n_c(t), # pending jobs n_p(t)
Explanation of the effect of new workers available after completion should be rephrased in terms of the T_compl_av, the average completion time.
Assume that s_max has upper limit s_lim. Then each job is finished after a finite time. How does this change long-time behaviour ?
What happens when the total number of workers w(t) is demand-regulated on time scales much longer than the longest site completion period ? d/dt w = (n_p-n_pGoal)/tau
A slightly relavant paper is Wilhelm03.

Paper-Project 2: PL fluctuations of chemical concentrations

2008-07-18T15:08:00.010+02:00

Max system:

Basic dynamics: d/dt y = + x y - x y^2 - c y
Fluctuations: x --> x + dx(t), multiplicative exp-correlated noise (tau_cor)
Inverted PL for P(x) in the limit of small tau_cor
Idea: P(x) corresponds to the stationary, long-time case. Then x(t) is effectively white. Small tau_cor means just small effective variance of the white x(t).
Check if Fokker-Plack-Eq. is solvable analytically or numerically.

A similar system and result is discussed for gen expression in Nacher06:

Basic dynamics: d/dt y = + p - x y.
Fluctuations: p --> p + dp(t), additive white noise. Seems to be not so important.
x --> x + dx(t), multiplicative white noise
Stationary limit of Fokker-Plack-Eq. is solvable for large conc. x and gives a PL.
PL is of normal decaying type

Fractional Fouriertransform

2008-07-18T11:08:00.001+02:00

The Fractional Fouriertransform can be viewed a rotation operation on the time frequency distribution.

Zanette's Nature-Essay "Playing by numbers"

2008-07-18T11:06:00.001+02:00

Nature 453, 988-989 (19 June 2008)

In 1955, social scientist Herbert Simon pointed out that Zipf’s law can be quantitatively explained by assuming that the usage frequency of a word increases proportionally to its previous appearances.
As the message flows, a context emerges, favouring the appearance of some elements at the expense of others.
Segmentation, for instance, can be used to detect portions of a sequence that differ as much as possible in the frequencies of different symbols. The product is a dissection of the sequence into domains which are maximally divergent.

Es existieren weitere Artikel zu ScienceAndMusic.

Goal of CSK dynamics

2008-06-24T17:13:00.008+02:00

We usually treat the spontaneous motion of CSK-bound beads as a random process. Could it also be viewed as goal-oriented behaviour ? Max Sajitz-Hermstein recently pointed out that the motion of whole cells with fractional superdiffusive MSD is interpreted by some researchers as an optimum foraging strategy.

On the other hand, the cell "thinks" of the bead as a substrate and excerts (power-law-correlated) forces onto this substrate, maybe in order to move itself along the opposite direction ?

It would be interesting to figure out the motion of a whole cell that results from the single stress fiber behaviour, as it was deduced from the bead diffusion measurements. We have to consider the effect of many stress fibers, probably acting independently from each other. Also the geometry is different. But with a broad distribution of the durations of the force pulses, there should indeed result a superdiffusive cell motion.

Weak Bead-CSK-Links

2008-06-18T17:17:00.011+02:00

Assume there is a preexisting CSK stressfiber network with a particular node (position r_N(t)). The node is performing a random motion with a fractional, superdiffusive MSD. A bead (position r_B(t)) is floating in a viscous medium (friction coefficient gamma), where it experiences white noise thermal forces (amplitude vN of velocity fluctuations), causing a diffusive motion. Now a weak spring (constant k) is connecting the bead to the CSK node. The visco-elastic relaxation time is tau=gamma/k. How does the MSD of the bead linked elastically to the node look like and how does it depend on vN and tau ?

The total force F(t) acting on the bead is given by:

At any moment, the forces from the node, from friction, and from thermal forces are balanced. Setting F(t)=0 yields a dynamic equation for the bead motion, which was solved numerically (program boreas_1).

Parameters:
NMX=64*1024;
dT=0.01;
plExp=1.5;
difus=5.0;

Here the MSD for small thermal noise:

And for larger thermal noise:

Note that depending on the binding strength (~1/tau), the combination of the diffusive, plateau-like and fractional superdiffusive powerlaws can produce an apparent sub-diffusive regime at intermediate lag-times.

The following figure shows the individual contributions to the total MSD:

New Lectures

2008-06-16T17:41:00.003+02:00

Complex Systems:

Biomechanics:

Dynamic crosslinkers (exercises, solutions)
Minimal sliding filament model (exercises, solutions)

Apollon-Project: Bead and Fiber SWD

2008-06-12T12:12:00.024+02:00

How is the SWD (form of distribution, mean value Dm, kurtosis kur) of the bead related to the SWD of the individual fibers ?

Assumptions:

Spider geometry with N_fib fibers.
Individual fibers have identical exponential SWDs with average step width dMean.

Program: Apollon_1

Parameters: dMean=1.0

Results:

dMean acts as a scaling parameter.
NFib acts as a shape parameter.
sqrt(NFib) ~ Dm ~ dMean
1/Nfib ~ kur <--indep--> dMean

As for the SWD itself, one can either analyze the distributions P(x), P(y) of the cartesian components of the bead displacement, or the radial distribution P(r).

In the case of P(x), one finds as a function of N_fib a gradual evolution from exponential to Gaussian:

The radial distribution P(r) corresponds, of course, to xP(x) and approaches in the limit of infinitely many fibers a Rayleigh distribution:

So far, the well stiffness kWell was a constant. What happens to the above results when kWell grows in proportion to the number of myosins NMyo in the system ? NMyo is proportional to the number of fibers NFib and to the average fiber size fMean, which depends in turn on the remodelling rate |df/dt| (fibers grow to higher peak sizes within the same time).

With myosin-depending stiffness, KWell ~ NMyo ~ NFib. Since dMean ~ 1/kWell, we get an extra factor dMean ~ 1/NFib. Together with the factor dMean ~ sqrt(NFib) we arrive at

dMean ~ 1/sqrt(NFib).

Note that dMean determines the overall diffusion constant D of the MSD. So with myosin-depending stiffness, D will decrease with the number of fibers. With constant stiffness it will increase. The experiments favor myosin-dependent stiffness.

Note also that kWell does not change the shape of the SWD, only the scale.

It is possible to reasonably fit the SWD to experimental data, at least within the main peak regime:

There is a certain flexibility: a higher N_fib can be approximately compensated by a smaller dMean.

[MT 9] Wiener–Khinchine theorem

2008-05-28T13:52:00.013+02:00

Reminder: We consider a real-valued, stationary random signal f(t) with zero mean (the average has been substracted, in order to yield the pure fluctuations).

The autocorrelation function (ACF) of f(t) is defined as:

The power spectral density (PSD, power spectrum) of f(t) is defined as:

where f_T(t) is the function f(t) restricted to the intervall [-T/2..+T/2] (i.e. all values outside the intervall are set to zero) and f_T(\omega) is the Fourier trafo of f_T(t).

According to the W.K. theorem, the ACF and the PSD form a Fourier pair:

The following properties are obvious:

The autocorrelation function C(\tau) of f(t) must be real.

C(\tau) should be an even function of lag time: C(-\tau)=C(\tau), because stationary signals should not behave different when we go forward or backward in time.

The Fourier transform of an even real function is, again, an even real function, so S(\omega) is even and real. The reality is clear, anyway, since S_f(\omega) is the modulus squared of f_T(\omega).

The value C(\tau=0) gives the statistical variance of f(t).

The value S(\omega=0) is proportional to the mean value of f(t).

Two cases are or special importance:

Exponential correlation:

Here, the variance is C_0.

For powerlaw correlations

where \tau is a dimensionless time and \gamma is between 0 and 1 one obtains:

Here, the variance C(\tau=0) and S(\omega=0) are divergent.

See for reference this and that link. These two references also contains interesting hints about the strange things that happen when \gamma approaches or even exceeds 1.

[MT 8] Comments to [MT 7]

2008-05-20T14:00:00.010+02:00

Assume again we have a Poisson train of assembly steps (AS) with average number k_av=10 per measurement time intervall. Let each FU require m=2 of the AS. Mentally, we can definde a super-Poisson process, where each pair of AS is merged to a FU-super-unit. Then one might naively think that

P(12,10)=P(having 12 AS done when there would normally be only 10)

might be the same as

P(6,5)=P(having 6 FU done when there would normally be only 5).

However, the Poisson distribution has no such scaling property:

P(12,10)=0.09478033009, but
P( 6, 5)=0.1462228081.

Thus, in general, P( k / m, k_av / m ) is not equal to P( k, k_av ), and the super-Poisson distribution P(n | m, k_av = R dt) defined in [MT 7] is not just a rescaled Poisson distribution.

Note also that the super-Poisson distribution has one extra parameter, m, which brings more flexibility in fitting the experimental SWDs.

[MT 7] Stochastic Assembly Processes

2008-05-16T14:13:00.014+02:00

In our biophysical model of spontaneous bead diffusion, the persistent motion of the bead is eventually traced back to selfassembly processes in the stress fibers. It is therefore important to consider the statistics of such processes in a simple model.

Let's assume that the goal of a biochemical factory is to produce Functional Units (FU) of some kind. The fabrication of each new FU requires a fixed number of m Assembly Steps (AS) to be performed (successively). The ASs are accomplished at an average rate R, with exponentially distributed inter-event-times (Poisson process). We are interested in the probability P(n | m, k_av = R dt) that in a given time intervall dt, n complete FUs are fabricated.

Note that the case m=1 is described by the Poisson distribution. The cases m>1 can be viewed as a generalization of the Poisson distribution. The intervalls between successive FUs are completed are described by the Gamma distribution.

It is straight forward to compute P(n | m, k_av = R dt) in a Monte-Carlo simulation, where exponentially distributed random durations tau_i are added, until sum(i=1..n) tau_i >= dt. Always when this happens, the resulting random integer n is recorded in a histogram and the counting process restarts. It is also possible to treat n as a continuous random variable by adding to the integer n the fractional part of the next FU that is already fabricated at the end of intervall dt.

[AD 18] sd-model: parameter dependence of SWD

2008-05-15T16:29:00.004+02:00

Starting from the fit parameters of the low b group, the effects of varying the three relevant parameters on the SWD was explored:

[AD 17] sd-model: fitting MSD & PSD

2008-05-15T13:48:00.011+02:00

Using the program dmsd.cpp and Carina's averaged data for the two groups of trajectories with low and high b values, a fit was attempted. The kurtosis was not considered.

Results for the low b group:

sigN2=1.2e-4;
beta=1.5;
lambda=0.05;
w1=145.0e-3;

Results for the high b group:

sigN2=0.6e-4;
beta=1.8;
lambda=0.35;
w1=20.0e-3;

With the above fit parameters, the simulated step width distribtions, evaluated for the smallest time intervall, look for the two groups as follows:

The SWDs are affected by the parameters lambda and w1 (for the continuous Poisson process) and by sigN2 (for the Gaussian background noise). The tails seem to be dominated by the Poisson process and are close to exponential.

[GR 2] Quantifying hopping-and-stalling motion

2008-05-13T08:44:00.018+02:00

Intuitively, hopping-and-stalling-motion (in the case of a discrete time random walk) is marked by

long chains of successive steps in which the particle remains bound within a box of small volume,
separated by a single long jump, carrying the particle to a new box.

The characteristic dimensionless numbers seem to be

The average number of steps within the small boxes
The ratio r=L/s of the long jump length L to size s of the small box

The problem is how to identify the long jumps without introducing an artificial threshold.

Let us first assume that the critical ratio r is prescribed. We can then start at the beginning of the trajectrory (t=0), go through all steps (t->t+1), and keep track of the box size s(t) covered so far (In 1D, the box size s(t)=x_max-x_min is the distance between the most right point x_max and the most left point x_min visited since the last jump). As soon as a single step exceeds s(t) by a factor of r, we have a jump. Then the N_b is known for this special chain and we continue with the next one. When we reach the end of the trajectory, we can compute the average N_b, L and s.

These averages, however, depend on the prescribed parameter r. If r is small, even not so long steps count as jumps and the typical chains will be short. Consequently, the averages of N_b as well as s should become small, and the same holds for L. Reversely, all averages increase with r.

Is there a way to define the "natural value of r", which is inherent to the trajectory ? If the above averages, when plotted against r, would for example show a sudden change of slope at some r_0, this could serve as a characteristic value, similar to the so-called "ellbow criterion" in cluster analysis.

It is also interesting to consider the special case r=1, where a jump must at least exceed the box size (jumping just out of the box): Starting at the beginning of a new chain, s(t) is increasing or staying constant at every step. Once it reaches a certain size, it becomes difficult for any single step to exceed s(t). Thus, (in particular when many successive steps are heading into the same direction, leading to a steady increase of s(t)) the r=1 criterion might never be met. In such a case there is definitely no hopping-and-stalling motion !

Note that for any trajectory of finite length, there is a maximum jump ratio r_max which is barely met a single time, separating the trajectory into two stalling clusters. We can define the hopping-and-stalling criterion to be fulfilled when r_max>1.

Unfortunately, it is not clear how to extend this idea to the 2D case.

[GR 1] Formalizing qualitative properties

2008-05-10T17:40:00.007+02:00

Sometimes the human mind intuitively recognizes a new property P of some object X, which can initially be stated in vague, qualitative terms only. For example: "This trajectory is marked by subsequent periods of hopping and stalling".

One would then like to define a mathematical function p(X) which returns for each suitable object X a number p that measures to which extend this property is really present - the "p-ness" of X.

In order for this to work, the objects X must be parametrizable by a set of N numbers, so that each concrete sample object occupies a point in an N-dimensional space. The function p(X) is a scalar field over that space.
The p-function should reflect our intuitive property as closely as possible.
By comparing p(X) with p(Y), one can order X and Y with respect to the p-property.
Any monotoneous function q(p) would also allow for the ordering, so there is some ambiguity in the definition.
If the values of p are known for some "standard objects", one can even make qualitative statements: "The p-temperature is larger than that of freezing water, but smaller than that of boiling water".
If the property is "binary" (as in the case of "subdiffusive/superdiffusive") one must provide a standard object defining the border line ("just diffusive").

[MT 6] Continuous Poisson distribution

2008-04-30T16:49:00.018+02:00

It would be useful to have a generalization of the Poisson distribution for continuous event numbers k. This can be achieved by replacing the factorial by the gamma-function:

The following plot compares the discrete PDF (black) with the continuous one (orange) for lambda=3:

The same for lambda=33:

Theoretically, the kurtosis of the discrete Poisson distribution is 1/lambda. The following plot compares the numerical kurtosis=f(lambda) in the discrete (orange) and in the continuous version (green) with the analytical expectation (black):

The continuous Poisson distribution is obviously related to the gamma distribution with a=k+1. However, in the gamma distribution with shape parameter a and rate parameter b,

it is x=lambda which is considered as the random variable.

[AD 16] Effect of microstep rate on MSD and PSD

2008-04-30T10:28:00.008+02:00

If the Poisson statistics is combined with the s*d process, the VAC depends on the system parameters as follows:

While w_1 merely affects the scale of correlated bead fluctuations, the number of microsteps lambda per discretization intervall enters non-linearly into the formula.

The following plot shows the effect of lambda on the MSD:

And, even more pronounced, on the PSD:

The parameters were as follows:

sigN2=0.5e-4; // variance of white noise [um^2]
beta=1.5; // PL exponent of MSD [1]
lambda=2.0; // aver. nr. of microsteps per dT [1]
w1=5.0e-3; // width of single microstep [um]

It is inctructive to compare this with the effect of changing w1 at fixed lambda. The MSD:

The PSD:

Obviously, the effects are rather similar. This comes in handy, when lambda is adjusted to fit the experimental kurtosis. Then the simultaneous changes in the MSD/PSD can be compensated via w_1.

[AD 15] Poisson-distributed SWD

2008-04-30T09:52:00.009+02:00

We now assume that the bead displacement d_n during time interval n is the sum of k_n unresolved microsteps, each of (average) width w_1:

Let the k_n fluctuate according to a Poisson-distribution with average event number lambda (given by some event rate times the discretization time interval):

which has a variance of

Then one obtains for the mean bead displacement per discretization time interval

and for the variance

[AD 14] Test of MSD and PSD in s*d model

2008-04-29T14:02:00.014+02:00

Within the sd-model, it can be shown that the VAC depends only on the mean value d_mean and on the variance d_var=sigD2 of the step width distribution, not on the detailed distribution function. The PL-exponent is not at all affected by the step widths.

As a test, the correlation function Css(m) of the sign factors was chosen as

This form assures a smooth transition to the required value of Css(0)=sigS2=1/4. The PL regime starts at m>>1. Then, the VAC of the correlated part of the random process yields

For simplicity, it was assumed for the evaluation that

Using the discrete MSD formula from [MT 6], one finds the expected behaviour:

A fast Fourier transform of Cvv yields the corresponding PSD (the frequencies are in Hz):

The simulation parameters were:

dT=0.1; // simulation time intervall
NMX=16*1024; // nr. of sim. time steps

sigN2=0.5e-4; // variance of white noise
dMean=0.5e-2; // mean step width
sigD2=dMean*dMean; // variance of step width
beta=1.25, 1.5, 1.75; // PL exponent of MSD

Note that a mean step width of about 5 nm seems to be a reasonable value, in order to get the sub- to superdiffusive transition time correct.

[MT 6] Discrete relation between VAC and MSD

2008-04-28T10:46:00.009+02:00

There are some subtle points to be considered when the continuous VAC-MSD formula is discretized. The following relation, however, should work:

A good test is to insert the VAC of differential white noise

into the above formula. One obtains for n>=1 the following result:

\!0)=2\sigma^2" alt="\mbox{MSD}_x(n\!>\!0)=2\sigma^2" border="0"/>
which is the correct flat plateau at a level proportional to the variance of the spatial white noise.

[PP 2] Start of a new paper manuscript

2008-04-24T13:55:00.013+02:00

Titel:

Cytoskeletal fluctuations as a stochastic process with four parameters

Abstract:

Information about biochemical processes in the cytoskeleton of living cells can be obtained by attaching microbeads to the cytoskeleton and observing their random motion. Sampling the two-dimensional trajectory in equidistant time intervalls generates a discrete time series (X_n,Y_n). From these raw data, a variety of statistical quantities can be computed, in particular the mean squared displacement (MSD), the velocity autocorrelation (VAC), the power spectral density (PSD), the step width distribution (SWD) and the turning angle distribution (TAD).

We show that all the above spectra and distribution functions can be consistently reproduced by a very simple stochastic process. It consists of a persistent random walk with additive white gaussian noise. The random walk has steps of finite variance, but with long-time (powerlaw) correlated directions. The total process is characterized by only four parameters, which are easily extracted from the data. While under stationary conditions the MSD, VAC and PSD are directly related to each other, it is demonstrated that the SWD is an independent property of the process. The experimentally observed leptocurtic SWD at small lag times is compatible with a Poisson distribution.

Based on the abstract stochastic process, we present a biophysical model of an evolving acto-myosin stress fiber network attached to the bead. Persistent phases of growth or degeneration of remodelling fibers give rise to a superdiffusive bead motion on longer time scales. Within each persistent phase, the discrete adding and removing of fiber building blocks leads to force steps that occur with Poisson statistics and thereby create a non-gaussian SWD. The characterization of bead trajectories by only four key parameters in connection with our biophysical model will be usefull to investigate the effect of pharmacological treatments of the cytoskeleton in the future.

To access the latest version of the manuscript, see the right sidebar and click
Manuscripts/4ParameterProcess

[AD 13] Interpretation of the force steps

2008-04-23T16:17:00.009+02:00

When observing the spontaneous bead fluctuations, we find that the trajectory consists of phases in which the bead moves persistently in one direction. During such a persistent phase, the kurtosis is positive.

Our physical interpretation is as follows: Since the bead is bound predominantly elastic, the persistent motion means that the force causing this motion is increasing more and more: a force ramp. However, the positive kurtosis indicates that the force increases not continuously, but in steps.

Note that what we need are not force pulses (going up and down), but steps (going up and remaining at the high level, until the next step increases the force even higher). Symetrically, a step-wise decrease of the force would move the bead in the opposite direction.

Biologically, the force steps should have something to do with the assembly/disassembly of acto-myosin stress fibers. My interpretation is as follows:

Upward/downward force steps occur each time when a new (force-generating) building block has been added/removed to one of the acto-myosin fibers.

This building block could be a bunch of myosin motors which have been added/removed to/from the pool of other motors that already contributed to the average fiber force before. It could also be a whole pre-assembled small acto-myosin filament that attaches itself in parallel to an excisting filament bundle.

It is important that the motors are working, at least collectively, in a "processive mode": There should at any moment be a finite average force level and never should all motors simultaneously let go the actin filament.

If we take the "fit parameters" from [AD 12] seriously, in order to achieve in one step a bead displacement of 10 nm, assuming a elastic binding stiffness of 1 nN/um, we need a force step of 10 pN. This is more than what a single myosin motor can probably achieve.

[AD 12] Results for the SD model

2008-04-23T15:40:00.006+02:00

The SD model in two (independent) spatial dimensions has been implemented and it seems to work.

According to the SD model (see [AD 11]), the autocorrelation should have a powerlaw distribution with some prescribed exponent p. This will produce long persistent chains of the same sign/direction. The length distribution of the chains will be a PL with another exponent gamma. The analytical relation between p and gamma is not known yet. In the simulation, it is actually the exponent gamma which is prescribed.

For the step length d, a Poisson distribution was chosen. Parameters are the width wid1 of a single step and the average time tauAv between successive steps.

Another parameter is the variance sig2Noise of the additive white Gaussian noise.

For the choice:

wid1=1.0e-2;
tauAv=1.5;
gamma=-2.5;
sig2Noise=0.25e-4;

and a maximum chain length of 1e4, the MSD looks as expected:

Also the kurtosis behaves as in the experiments:

The step width distributions (without noise) show the characteristic exponential-like form for small lag times:

It is interesting that the "fit" parameters indicate force steps with an individual magnitude of 10 nm, following in average intervalls of 1.5 sec. This should be intrepreted biologically.

[PP 1] Lab Meeting Presentation April

2008-04-22T08:01:00.011+02:00

Analytic theory of anomaleous diffusion:

The kurtosis thriller:

"Realistic" models:

ZigZag model

Nonlinear Rheology:

[AD 11] Solution of the PL-Kurtosis puzzle - The SD model

2008-04-21T11:22:00.016+02:00

The question if powerlaw diffusion is compatible with positive kurtosis and, in particular, with an exponential step width distribution, could finally be settled by directly constructing a discrete time series with the requested properties:

We consider a simple stochastic process, which produces a time series of values x_k=x(t_k) for discrete time points t_k = k*dt:

where the s_n are sign factors

with equal distribution function

and a powerlaw auto-correlation:

The d_n are positive step widths

with an arbitrary distribution function of finite variance (Gaussian, Exponential, Poisson...)

but without any temporal correlations (white noise):

which means that

We assume there are no cross-correlations between the sign-factors and the step widths:

The velocity of the nth step is

For the auto-correlation function of the velocity we obtain

which can be factorized because the s and d and uncorrelated:

which yields for all m except 0:

This means that the process has powerlaw VAC (and PSD and MSD) for arbitrary step widths distributions.

By using separate processes x_n and y_n of the above kind for the different spatial directions, the method can be generalized to multiple dimensions.

[NR 2] Chain of NKBs

2008-04-19T11:09:00.010+02:00

We first consider a chain of N linear Kelvin bodies (LKBs):

Assume that all dashpots are equal, but that the spring constants k_n differ for each body n.

Since all LKBs feel the same force F(t), a separate evolution equation can be written for each individual body:

We are interested in the total length X(t) of the chain, which is given by:

The total microstate of the system is described by the vector of the N individual lengths x_n(t), but externally only the global variable X(t) is relevant.

In this linear case, we can define a response function S_n(w) for each individual LKB:

Remarkably, we can even define a global response function of the LKB chain in frequency space, which directly connects the total chain length X with the external force F:

This means that in the linear case, we don't have to keep track of the N microvariables if we just want to know the global X<-->F response.

Now we replace the LKBs by the nonlinear Kelvin bodies described in [NR 1]. All we can do is to write down the separate DGLs for each individual NKB:

There is no abbreviation in this nonlinear case to obtain X(t) from F(t): We have to solve the N coupled differential equations simultaneously and at each time sum over the individual lengths.

So this simple example demonstrates:

We cannot in general treat a non-linear system like a black box and write down response equations which depend only on the externally accessible macro-variables. We rather have to use an evolution equation for the full multi-component micro-state, which will boil down to the solution of a coupled set of nonlinear differential equations.

[NR 1] Nonlinear Kelvin Body (NKB)

2008-04-19T10:22:00.013+02:00

Consider a Kelvin body, i.e. a spring in parallel with a dashpot:

The body is externally described by its length x and force F.

Assume that the rest length (at zero force) is zero. The dashpot has the property:

A linear spring would obey:

The two forces add up to the total force of the Kelvin body:

We can solve for the temporal change of x to obtain a diff.eq. describing the evolution of x(t) from its starting value for any given external force F(t):

A Fourier transformation yields

which allows to define a linear response function S(w) of the system:

Let the spring now have some nonlinear force-length-relation N(x):

The evolution equation then reads:

This nonlinear diff.eq. can also be solved (at least numerically) for any applied F(t). However, it is in general not possible to define a response function. This will have important consequences when we consider more complex systems build from many Kelvin bodies.

[AD 10] Can pure PL vel. correlations produce exp. step width distributions ?

2008-04-17T17:45:00.012+02:00

We have seen before that pure PL velocity correlations can produce positive kurtosis, however so far only in a very unrealistic way. Truely convincing would be a PL trajectory with an exponential (in stead of Gaussian) step width distribution (SWD). Can the velocity phases be tuned in such a way ?

To answer that question, en ensemble of 1000 trajectories was generated, with velocity phases randomly distributed between -pi an +pi. For each trajectory, the SWD was computed and least-square-fitted to an exponential distribution. The resulting error was high for all trajectories and fluctuated between about 20 and 40 (a good fit would have an error smaller than 1).

The trajectory with the lowest error was nr. 558. However, even this guy had still an extremely Gaussian SWD:

The kurtosis-versus-lagtime was also computed for 558 and compared to that of an artificial trajectory with truely exponential steps:

Obviosly, 558 is still far away from our goal.

In conclusion, it seems very probable that pure PL trajectories cannot have exponential SWDs. Or is the island of exponential-SWD concentrated in a very tiny region of the huge phase space, so that even among 1000 random probes not a single one comes close ?

Of course, PL trajectories might be compatible with other types of realistic distributions which positive kurtosis, like Gaussian peaks with abnormally fat tails.

[MT 5] PSD of ramp-like random processes

2008-04-17T10:08:00.003+02:00

In the same way as [MT 4] we can consider single "steps" of a different form. Most relevant for our CSK fluctuation models are sawtooth-like pulses of the form

2&space;&space;&space;\end{cases}" alt="s(t)=\begin{cases}
0 & \text{ if } x<0 \\
t & \text{ if } 0\leq x\leq 1 \\
2-t & \text{ if } 1\leq x\leq2 \\
0 & \text{ if } x>2
\end{cases}" border="0"/>

The Fourier transform reads

I didn't perform the calulations analogous to [MT 4] yet, but probably the asymptotic power spectrum would now be like 1/w^4 for large frequencies (short times).

[MT 4] PSD of step-like random processes

2008-04-16T19:57:00.043+02:00

In this post we construct a rather general "step-like" random process, consisting of many superposed steps of widely differing lengths tau_i and heights h_i. We are interested in the ensemble averaged power spectral density (PSD) of the process.

First we define a single step (more accurately: box) function of length tau and height 1:

.

The Fourier transform of b(t) is

Next we define a Poisson spike train with average spike rate lambda:

This process has a frequency independent PSD, equal to its rate:

We next generate a superposition of infinitely many positive, temporaly shifted, equally long and high boxes by convoluting the single box with the spike train:

Note that when the inter-box-intervall 1/lambda is much longer than the single box duration tau, successive boxes will typically not overlap.

The Fourier transform of random process f(t) is

and for the PSD we obtain

which yields in our special case

Our present random process contains only positive boxes. For a more symmetric process, we use two independent Poisson spike trains with the same average rate, but with opposite signs:

Since the two spike trains are statistically independent, the PSD of the term in brackets in just the sum of the two individual PSDs:

Thus, the PSD of the symmetric process is

.

The above process consists only of boxes of the same duration tau. We next consider a superposition of many uncorrelated sub-processes with different tau_i. Again, the total PSD will simply be the sum of all sub-processes. Let the density distribution of the tau_i be Q(tau). We assume further that the range of tau_i is restricted

and that Q(tau) is a "smoothly varying function" within the above limits. The total PSD is a weighted integral over all sub-processes:

One therefore obtains

.

For very small frequencies,

the sin^2 function can be replaced by the square of its argument,

yielding

where the remaining integral gives just a constant factor c_1.

For very high frequencies,

the sin^2 is a rapidly oscillating function, multiplied by a slowly varying Q(tau)-function. It can then be replaced by its mean, averaged over one oscillation period:

We thus get

where the remaining integral gives another constant factor c_2.

In conclusion, the PSD of step-like random processes with a limited range of step durations is constant (white noise like) for small frequencies (long times) and decays like 1/w^2 for large frequencies (short times).

Accordingly, if the above step-like process would correspond to the force fluctuations within an elastic medium, the bead MSD would have an asymptotic plateau for times much longer than the upper time constant tau_2. Unfortunately, the interesting w^-2 part of the force fluctuations at short times is not easily observable due to the noise floor.

The model could easily be extended to include steps of varying height as well. This would, however, not change the asymptotics of the power spectrum.

The above calculation does, of course, not prove that all processes with a w^-2 power spectrum are step-like, but at least the comment in the Bursac paper appears more reasonable now.

[MT 3] Differential Discrete White Noise

2008-04-12T10:35:00.037+02:00

Assume a white noise discrete time series (sampled at intervalls of dt) with variance sigma^2, here interpreted as a fluctuating position variable:

The discrete velocity is defined as:

Its autocorrelation yields (using the Cxx(m) above):

The PSD of the velocity fluctuations is given by the Fourier transform of Cvv. We have to use the Discrete Fourier Transform (DFT) here:

The DFT is based on discrete times

and the corresponding frequencies

where the indices run as follows:

One obtains the following PSD of the velocity fluctuations:

In a log-log plot, this looks as shown in this figure. The slope is always 2. So in the case of the experimental bead PSD, the apparent fractional slope of the high frequency wing comes from the superposition with the low frequency fractional powerlaw.

For a detailed derivation of the above formulae, see
SharedDocuments/Notes/DifferentialDiscreteWhiteNoise.pdf

The following figure shows a comparison between the above theoretical PSD and experimental data, measured for a purely noisy bead with a flat MSD:

[MP 1] Version "MultiLevel6"

2008-04-11T15:49:00.011+02:00

Collective Relaxation

In each relaxation step, all nodes are now moving SIMULTANEOUSLY along there negative gradient vectors.
The step width of each node is proportional to the magnitude of its gradient, i.e. nodes move different amounts.
In my experience, this strategy avoids problems connected with the wrong order of successive node movements.
The numerical effort of the collective optimization is the same as with successive optimization.

Naive Optimization Strategy

The global gradient ( = derivatives of the mismatch potential with respect to all node coordinates) is computed.
All nodes move simultaneously a small step.
The global gradient is computed again and so on.
Normally, line search is done in a different way: First walk in one direction until no more improvement is possible, and only then choose a new direction.
For some reason, the naive approach seems to work better. Don't really understand why.

New Graphics Function

The original intensity distribution is plotted in red.
The momentary fit is plotted in blue.
The goal intensity distribution is plotted in green.
All colors are plotted on top of each other, with slightly different pixel sizes.
So, in the end, the blue fit should coincede with the green goal.
Graphics windows can be advanced by clicking with the mouse somewhere inside.

Performance of Multilevel-Version-6

Amazingly, many test problems created with GenerateRndImages() or GenerateFullRndImages() (see [7]) are solved nicely, without any elastic forces.
However, the amount by which the mismatch potential can be relaxed, depends on the individual case.
Here, also the exact setting of the step width in the relaxation routine has an influence.
The real cell images with LoadRealImages() are not working at all (- so far..) .
Elastic forces don't seem to improve the results.

Elastic Forces

Are implemented.
Pure elastic relaxation on single grid level has been tested.
Combination of mismatch and elastic forces don't really seem to work well yet.

Random Test Images

GenerateRndImages(128,100,0.0,1.5) : Generates an orig and goal picture with resolution 128x128 pixels. Pixels are dark by default, but 100 random pixels have a finite random brightness. In the goal image, those pixels (more precisely: the underlying nodes) are shifted into a random direction by a random amount between 0.0 and 1.5 pixel sizes.

GenerateFullRndImages(32,0.0,1.5): Similar to above, but here ALL 32x32 pixels will have some random brightnesses.

To generate a new random case, change the integer variable SIT in DefParameters().

Idea: Separate Registration and Traction Computation

If the mismatch potential can really be relaxed without any elastic forces, this would be numerically effecient for images with many dark pixels.

Then, after we know where each bright node has moved to, the elastic problem and force calculation could be done as a second step.

Simplified Point Spread Function

Since the new minimization strategy does not need second derivatives, it was possible to go back to the lower-order, triangle-like PSF, leading to a second order B-function.

COI based Coarsification and Refinement

Lolo's idea is used, i.e. in the coarsification step each low-level supernode is placed at the center of intensity (COI) of its associated high-level component nodes.

In the refinement step, the group of component nodes is shifted, as a whole, by the same amount their supernode has been shifted in the grid below.

Coarsification, low-level relaxation and refinement together now comprise a true correction: Low-level and high-level information are combined.

It is ensured by the COI method that mere rigid shifts between the orig and goal images are treated exactly.

With the COI method, low-level grids are not uniform any longer. Thus, elastic springs have different rest lengths.

It is not clear yet how to implement a given Poisson ratio (requiring different spring constants for vertical/horizontal and diagonal springs) on such a non-uniform grid.

Beyond Powers of 2

Ben had the idea to use more flexible schemes for the multigrid method coarsification, so that the linear picture sizes doent need to be powers of two each time.

For example, one could also use 3x3 groupings instead of 2x2.

Or, one could be satisfied with a lowest grid resolution of, say, 5x5, instead of 1x1, but still use the normal 2x2 refinement.

[AD 9] "ZigZag" - Model of CSK remodelling

2008-04-09T10:07:00.025+02:00

General concepts:

We are interested in the spontaneous motion of CSK-bound micro-beads.

In micro-particle rheology, it is common to explain the motion of the test particle by the interplay of two factors: First, there are external forces acting on the particle. Second, there is a rheological response function, which determines how the particle moves for a given external force profile.

In our case of the cell, this distinction is not so obvious, since the bead is bound to an acto-myosin fiber network with active and passive components. The same network plays the role of a medium and also provides "external" forces, due to remodelling events.

To make the conceptual problem clearer, let us imagine the fiber network as a mesh of linear springs with given spring constants and rest lengths. The bead is just a special node in this network. Its position is at any time at the mechanical equilibrium position.

Since living biological systems are due to constant turnover and remodelling, the spring constants and rest lengths are temporally changing. In particular, stress fibers can contract (shrinking rest length) and thus create prestress in the network. Consequently, the beads equilibrium position is constantly shifted around. This resembles the observed spontaneous bead motion.

Note that it is not clear in the above example how to separate the parts played by the medium and the forces. All we can say is that the parameters of our network are gradually changing and this is accompanied by shape changes of the network. It evolves through a succession of mechanical equilibrium configurations. After each remodelling "event", the prestress as well as the stiffness of certain fibers may have changed simultaneously.

However, to simplify the analysis, we nevertheless stick to the medium/forces distinction. In particular, we will treat the forces in more detail than the medium in the following model.

The properties of the medium can in principal be measured independently, namely by applying known forces (for example: a force ramp) to the test particle and observing the resulting motion. This must be done quickly enough to avoid remodelling during the test. Such experiments have been done and it was found that the CSK is predominantly elastic on such relatively small time scales.

For further simplification, we describe the medium by a radially symmetric potential well (reasonable in the sense of a statistical average). The bead would rest in the center of the well, if the prestress forces of the fibers would not change with time.

However, we should at least account for the correlated changes of stiffness and prestress in the remodelling network. As a reasonable approximation, one could set the radial stiffness coefficient of the elastic well at each moment to a value that is proportional to the momentary total prestress of the stress fiber network (sum over all individual fiber forces).

In reality, the network may consist not only of active, contractile fibers, but might also contain a passive sub-network. The stiffness contribution of the passive sub-network might be independent from the remodelling of the active stress fibers. In the following model, we describe the passive contribution by a constant, minimum stiffness of the well, which remains even if all active fibers disappear. This passive constribution is however chosen much smaller than the average active constribution in the steady state.

Besides the prestress forces of the fibers, there are also other, quite trivial reasons for random bead fluctuations in a real experiment, as for example noise in the detection system and thermal fluctuations. They can be described by white Gaussian noise of suitable variance, uncorrelated in the x- and y-directions. This positional noise is simply added to the simulated bead trajectories in the end.

Assume for moment that the well stiffness is constant. Then the bead position is at any moment proportional to the sum of all prestress forces acting on the node=bead. A bead moving with constant velocity in the direction of a particular fiber means that the force of this fiber increases linearly with time.

Since the MSD of the bead is a fractional powerlaw (at least above the noise floor), the fluctuations of the prestress of each fiber have to be long-time correlated, or must at least contain a very wide spectrum of short, medium, and very long lasting processes.

For simplicity, we assume that the prestress of each fiber grows with a fixed average rate for a certain duration. Then it decreases with another rate until zero. A new triangular cycle starts immediately, however with a new random duration, drawn from a wide probablity distribution.

Just as each cycle is statistially independent from the previous one, the triangular fluctuations of different fibers are also independent from each other. The fiber directions remain fixed in the model.

The model as described so far can explain many features observed in the experiments, but not the positive kurtosis of the bead's step width distribution at small lag times. A positive kurtosis is obtained when we assume that changes of a fiber's prestress do not proceed continuously, but in a "quantized" manner, i.e. in discrete units.

This is actually not unreasonable: A stress fiber can add a g-actin monomer to the end of a polymerizing filament. This increases the overlap between parallel filaments and offers new opportunities for myosin motors to crosslink. As a consequence, the prestress increases by a certain discrete amount (Note that the increase must not be prescisely the same each time).

The adding or removing of structual, force-generating units to a stress fiber is a stochastic process. It can be most easily modelled by a Poisson process: The regulatory system of the cell prescribes only the average rates of the assembly/disassembly processes. The actual time of each single event is random with exponentially distributed inter-event times.

The average assembly and disassembly rates are themselves changing, however in a more gradual manner. For simplicity, we assume that the disassembly rate is constant. The assembly rate is up-regulated to a value exceeding the disassembly rate during the growth phases of the fiber. It is down-regulated during the shrinking phases.

The model in detail (model parameters in red):

NFib fibers are attached to the bead in radial spider geometry with statistically isotrope directions.

The bead is bound in a radially symmetric elastic potential well with a prestress-dependent stiffness kWell(t).

The fibers exert a temporaly changing prestresses F_n(t)=f_n(t)*unitvec_n onto the bead (reflecting growth and shrinking of the fibers by biochemical processes).

The forces add up vectorially to F_sum(t).

The bead's mechanical equilibrium position R(t) is at each moment given by F_sum(t)/kWell(t).

Here, kWell(t) = kWell_Min + c * [sum_n f_n(t)].

The fiber growth times forces are log-normal distributed: tGrow_n = 1sec * 10^GaussRnd(eMean,eFluct).

Fibers grow with constant rate for a time tGrow_n. Then they decay with the same rate to zero force. Then another cycle starts with a new random tGrow_n.

The above temporal changes of fiber force f(t) are generated by an imbalance between a constant background disassembly rate lamdaDis (unit: Newton/sec) and a varying assembly rate lambdaAss(t).

In the model, lambdaAss(t) is switching between values 0 (shrinking phase) and lambdaAssMax (growing phase).

The fibers grow/shrink with Poisson kinetics: Fibers change size in discrete units. Each additional unit increases the fiber force by an amount frcStp.

A Gaussian white noise of variance sigma2Noise it added a posteriori to the bead position to simulate thermal and instrument noise.

Example of resulting MSD with and without (smooth) Poisson kinetics:

And the kurtosis:

The size of fiber 1 versus time:

Simulation-Parameters:

NMX=10; // nr. of fibers
lambdaDis=10.0; // disassembly rate in [pN/sec]
lambdaAssMax=2.0*lambdaDis; // on value of assembly (fiber force increase) rate
eMean=0.0; // distribution of fiber growth times ,
eFluct=1.5; // tGrow_n = 1sec * 10^GaussRnd(eMean,eFluct)
frcStp=1.0; // force increase per assembled unit in [pN]
kStiff=100.0; // elastic stiffness in [pN/um]
sigma2Noise=1.0e-4; // white noise force amplitude (plateau height) in [um2/sec]

[AD 8] Additive Gaussian Noise : Effect of Cor.Time

2008-04-09T08:57:00.007+02:00

In some of the more ballistic experimental data, the VAC seems to indicate an exponentially decaying correlation (probably followed by the expected PL decay, which might be on a too low level to be observable. The PL can be clearly seen in the PSD and MSD, however). Could this indicate that the additive noise is actually not white, but exponentially correlated with a characteristic time tau?

Simulated effect of varying tau on the MSD:

And the effect on the kurtosis:

Note that the used noise has a Gaussian PDF, so there should be no leptocurtic behaviour (The veloity phase distribution is a box [-pi,+pi].

[AD 7] Unrealistic Leptocurtic Mechanism

2008-04-07T19:09:00.009+02:00

Unfortunately, it turned out that the narrow Gaussian velocity phase distribution achieves the leptocurtic behaviour by an unrealistic mechanism: The momentary velocity (in x- and y-directions) has a strong peak exactly at the center of the simulation time interval:

This leads to a huge step in the displacements:

This in turn leads to some isolated contributions in the far tail of the step width distribution, which causes the positive cutosis.

This is not what we are looking for.

[AD 6] Gaussian Velocity Phase Distributions

2008-04-07T11:13:00.008+02:00

The result [AD 4] does not mean that PL-VAC cannot produce positive kurtosis. To check that, noise is completely turned off for the time being.

Starting from a particular trajectory with pure PL-VAC that showed a slightly positive kurtosis, the phases have been optimized in an evolutionary procedure (mutation + selection). By inspecting the increasingly leptocurtic trajectories, it turned out that their phase distribution evolved from the initial box [-pi,+pi] to something narrower and Gaussian-like.

The following graph shows kurtosis vs lag-time for Gaussian phase distributions, centered around zero and with varying FWHM:

That means that leptocurtic behaviour can be adjusted simply by the FWHM parameter. So our stochastic process is characterized by 4 parameters totally.

[MT 2] Relations between PL exponents

2008-04-07T10:08:00.005+02:00

Realistically, the effective medium of the cell should be described by a frequency dependent eleastic modulus

For simplicity we assume a purely elastic medium with x=1, so that the bead displacement is directly proportional to the momentary sum of forces. Then the PL exponents are related as follows:

Mean squared displacement:

Velocity Autocorrelation:

Velocity power spectral density:

Force power spectral density:

In our case we get for diffusive motion lambda=2 and for ballistic motion lambda=3. According to

P.Bursac et al., Nature Materials 4, 557 (2005)

a lambda of 2 is characteristic for force fluctuations that are finite but discontinuous (a series of force steps), whereas lambda=4 corresponds to uniformly continous fluctuations.

This could be checked by evolving the random phases in v(w) and then analyzing the resulting x(t) ~ F(t).

[AD 5] New model: "PL VAC + Gaussian Noise + X"

2008-04-06T10:52:00.005+02:00

Since we are not interested in shaped noise, the following concept would be more transparent:

The PL-VAC creates a trajectory with rich internal structure (a succession of more diffusive and more ballistic phases) and yields the correct MSD/VAC/PSD at medium and long times. It also makes the kurtosis slightly negative at long times.

Gaussian White Noise is added to the trajectory at the end. It only serves to produce a plateau in the MSD.

Another, yet unidentified Process X is reponsible for the positive kurtosis at small times.

Process X should have the following properties:

No significant distortion of the MSD/VAC/PSD
Produce increasingly non-Gaussian step width distributions with positive kurtosis as the lag time become small.
Candidate: Modulated Poisson Statistics.

[AD 4] Positive kurtosis comes from shaped noise

2008-04-06T10:16:00.006+02:00

We have shown before that a narrowed phase distribution yields a positive kurtosis at small lag times. It remains to be shown if the noise or the PL part is responsible for this behaviour.

Following graph: Kurtosis vs lag-time for wide phase distribution, when only noise or only PL is present:

Pure noise has no kurtosis and pure PL has a kurtosis decaying from zero to negative values.

Following graph: The same for narrowed phase distribution:

So it is clearly the noise part that creates positive kurtosis. The noise is shaped by the narrowed phase distribution in a way that makes the step width (= average velocity * time step) distribution non-Gaussian (compatible with the VAC).

This conclusion is reasonable: In the range of small lag times, the noise part anyway dominates over the PL part. Therefore, the PL-VAC cannot be responsible for the positive kurtosis.

[AD 3] Kurtosis in the "Noise+PL_Cvv" model

2008-04-04T16:30:00.003+02:00

When the phases of the velocity oscillations phi_k=arg(v(om_k)) are drawn independently from an equal distribution between 0 and 2 pi, the average kurtosis of the step width distributions seems to be zero (at least for small lag times, where the statistics is more reliable).

When the phase distribution is limited to a smaller intervall like [-0.8pi, +0.8pi], the kurtosis quickly raises to a high positive value at small lag times. From there it decays to slightly negative values at high lag times. This resembles the experimental data.

In conclusion, the full ensemble of trajectoties of the "Noise+PL_Cvv" model containes also a sub-ensemble which is fully compatible with the experimental data.

[MT 1] How to combine x- and y-step widths for the calculation of the kurtosis

2008-04-04T14:30:00.004+02:00

Situation:

We have a list of two-dimensional points (x_i, y_i), describing a particle trajectory.
This defines a list of steps (dx_i, dy_i)=(x_i-x_i-1 , y_i-y_i-1 ).
The steps have a certain 2D statistical distribution P(dx_i, dy_i).
We are interested in "The kurtosis of the step width distribution".

-----------------------------------------

Method 1: Merging dx- and dy-lists

Some people proceed as follows:

Show that the dx_i and dy_j are statistically independent.
Merge the dx_i and dy_i into a single list (ds_i) = (dx_1, dx_2,..., dx_N, dy_1, dy_2,..., dy_N).
Compute the kurtosis of the combined list (ds_i).

A simple example demonstrates that this procedure is incorrect in general:

Assume that P(dx) = Gaussian[mean=0,var=sigmaX].
Assume that P(dy) = Gaussian[mean=0,var=sigmaY].
Assume that sigmaX >> sigmaY, so that the 2D distribution is elliptically streched along the x-axis.
The distribution P(ds) of the combined list entries (ds_i) looks as follows: A narrow Gaussian peak riding on top of a broad one.
This pronounced central peak with broad wings results in a positive kurtosis.
But in reality the kurtosis of uncorrelated Gaussian random numbers is zero.

-----------------------------------------

Method 2: Separate kurtosis

The easiest solution is to give up the idea of a combined kurtosis and to determine two separate kurtosis for the dx- and dy-data. Finally the average may be taken - or not.

-----------------------------------------

Method 3: Using Euklidian distance

The width of step i is defined as

It appears natural to use the distribution P(dr) for computing "The step width kurtosis". However, in a 2D plane points at radius r around the center are geometrically over-represented by a factor of r,

and so a naive numerical histogram P(dr) will be distorted by this geometrical factor. The operation P(dr)-->P(dr)/dr fixes the problem in principle, but multiplies the noise at small dr. This is especially critical for computing the curtosis afterwards.

-----------------------------------------

[AD 2] Results for the "Noise+PL_Cvv" Model

2008-04-04T10:30:00.007+02:00

Here an example for three PSD curves, generated for different parameters b = 1.1 (black) / 1.5 (orange) / 1.9 (green) :

Here the corresponding MSD curves:

And here the trajectories:

[AD 1] The "Noise+PL_Cvv" Model

2008-04-04T09:30:00.007+02:00

Question:

Are Carina's data compatible with the simple assumption of strictly powerlaw-correlated velocity fluctuations on a white noise floor, or do we need extra model ingredients ?

Ansatz:

On a discrete time-/frequency grid, we generate a power spectral density (PSD) according to above assumption:

According to the Wiener-Kinchin theorem, P(w) is the Fourier transform of the velocity auto-correlation function (VAC):

The VAC in turn determines the mean squared displacement (MSD) of the particle:

In other words, the spectra PSD, VAC and MSD express the same information in different ways and are in the above model characterized by only three parameters a,b and c.

The square root of P(w) gives the frequency-dependent modulus of the velocity oscillations:

If we multiply |v(w)| by an arbitrary, frequency-dependent phase function, this will leave the PSD, VAC and MSD invariant:

But this operation will have a drastic effect on the shape of the resulting trajectory in real space and on other properties which are independent from the MSD:

In other words, we can explore a huge ensemble of different trajectories (by chosing for example random phases at each discrete frequency), which however all produce the very same MSD.

We can then find out if this space of trajectories is sufficiently rich to reproduce all the different signatures found in the data.