Simple pattern formation with cellular automata

A cellular automaton is a dynamical system where space, time and dynamic variable are all discrete. The system is thus composed of a lattice of cells (discrete space), each described by a state (discrete dynamic variable) which evolve into the next time step (discrete time) according to a dynamic rule.
\begin{equation}
x_i^{t+1} = f(x_i^t, \Omega_i^t, \xi)
\end{equation}
This rule generally depends on the state of the target cell $x_i^t$, the state of its neighbors $\Omega_i^t$, and a number of auxiliary external variables $\xi$. Since all these inputs are discrete, we can enumerate them and then define the dynamic rule by a transition table. The transition table maps each possible input to the next state for the cell. As an example consider the elementary 1D cellular automaton. In this case the neighborhood consists of only the 2 nearest neighbors $\Omega_i^t = \{x_{i-1}^t, x_{i+1}^t\}$ and no external variables.

In general, there are two types of neighborhoods, commonly classified as Moore or Von Neumann. A Moore neighborhood of radius $r$ corresponds to all cells within a hypercube of size $r$ centered at the current cell. In 2D we can write it as $\Omega_{ij}^t = \{x^t_{kl}:|i-k|\leq r \wedge |j-l|\leq r\}\setminus x^t_{ij}$. The Von Neumann neighborhood is more restrictive: only cells within a manhattan distance of $r$ belong to the neighborhood. In 2D we write $\Omega_{ij}^t = \{x^t_{kl}:|i-l|+|j-k| \leq r\}\setminus x^t_{ij}$.

Finally it is worth elucidating the concept of totalistic automata. In high dimensional spaces, the number of possible configurations of the neighborhood $\Omega$ can be quite large. As a simplification, we may consider instead as an input to the transition table the sum of all neighbors in a specific state $N_k = \sum_{x \in \Omega}\delta(x = k)$. If there are only 2 states, we need only consider $N_1$, since $N_0 = r – N_1$. For an arbitrary number $m$ of states, we will obviously need to consider $m-1$ such inputs to fully characterize the neighborhood. Even then, each input $N_k$ can take $r+1$ different values, which might be too much. In such cases we may consider only the case when $N_k$ is above some threshold. Then we can define as an input the boolean variable

\begin{equation}
P_{k,T}=\begin{cases}
1& \text{if $N_k \geq T$},\\
0& \text{if $N_k < T$}.
\end{cases}
\end{equation}

In the simulation you can find here, I considered a cellular automaton with the following properties: number of states $m=2$; moore neighborhood with radius $r=1$; lattice size $L_x \times L_y$; and 3 inputs for the transition table:

  • Current state $x_{ij}^t$
  • Neighborhood state $P_{1,T}$ with $T$ unspecified
  • One external input $\xi$\begin{equation}
    \xi_{ij}=\begin{cases}
    1& \text{if $i \geq L_x/2$},\\
    0& \text{if $i < L_x/2$}.
    \end{cases}
    \end{equation}
  • Initial condition $x_{ij} = 0 \; \forall_{ij}$

For these conditions a deterministic simulation of these conditions yields only a few steady states: homogeneous 1 or 0, half the lattice 1 and the other 0, and oscillation between a combination of the previous.

One possibility would be to add noise to the cellular automaton in order to provide more interesting dynamics. There are two ways to add noise to a cellular automaton:

The most straightforward way is to perform the following procedure at each time step:

  • Apply the deterministic dynamics to the whole lattice
  • For each lattice site $ij$, invert the state $x_{ij}$ with probability $p$

This procedure only works of course for $m=2$. In the case of more states there is no obvious way to generalize the procedure and we need to use a proper monte carlo method to get the dynamics.

A second way is to implement a probabilistic cellular automaton. In this case the transition table is generalized to a markov matrix: each input is now mapped not to a specific state but rather to a set of probabilities for a transition to each state ($m$ probabilities). Naturally for each input these sum to one. In this case we have $m$ times more parameters than before.

Solving stochastic differential equations with theano

In systems biology we often need to solve diffusion equations of the type $$df = f(x,t) dt + g(x,t)dW$$ where W is a white noise process; they’re the most common example of a stochastic differential equation (SDE). There are only very few cases for which we can analytically solve this equation, such as when either f or g are constant or just depend linearly on x. Of course most interesting cases involve complicated f and g functions, so we need to solve them numerically.

One way to solve this is to use a straightforward variant of the Euler method. You just need to do the usual time discretization $\delta t = T/N$ (with T the total time and N the number of steps) and then update your current value with the deterministic contribution $\delta t \times f$ and the stochastic contribution $\sqrt{\delta t} \times g$. If you are not familiar with stochastic calculus you may be wondering why is the time step multiplier $\sqrt{\delta t}$ for the stochastic part.

Continue reading “Solving stochastic differential equations with theano”

Simulating tissues with pressure

One small project I did was to code up a simulation of a growing tissue which feels pressure and where each cell has a dynamic state which depends on its neighbors and the pressure it feels. The idea is to reproduce some essential properties of morphogenesis. You can look at the code here. I am going to talk about the most interesting parts of the code.

I initialize the cells in an ordered lattice, with random perturbations in their positions, except those which are in the borders (bottom, left, right). Those are static and do not evolve in the simulation like the others. They are there just to represent the pressure from the rest of the body (huge) on the simulated tissue (tiny). This is not a very realistic assumption because many developmental systems have a size of the order of the body size, but we have to start somewhere!

All cells are connected with springs, which simulate adhesive and pressure forces in the tissue. If left alone, the system relaxes into a hexagonal configuration, since this minimises the spring potential energy. I integrate the harmonic oscillator equations using a fourth order Adams Moulton algorithm.

Now, it is important to realize that there are two time scales in the system: the pressure equilibration and cell lifetimes. We can assume the mechanical pressure equilibrates very fast, while cell divisions take their time. So what we do is run the oscillator system until equlibrium for each time step of the cellular state evolution, which we will talk about later.

Springs connect each cell. The larger colored circle is of the same size as the rest length of the springs. Thus, overlapping circles mean the spring wants to extend, while spaces mean the spring wants to contract. The color denotes the automaton state of the cell.
Springs connect each cell. The larger colored circle is of the same size as the rest length of the springs. Thus, overlapping circles mean the spring wants to extend, while spaces mean the spring wants to contract. The colors are explained below.

Continue reading “Simulating tissues with pressure”

Preprint review: Parameter Space Compression Underlies Emergent Theories and Predictive Models

So here’s a preprint I found really interesting [arxiv:1303.6738]. I’ll try to give a quick overview of the story in my own words.

The main concept used in the paper is the Fisher Information, which is no more than a measure of the curvature in the space of probability distributions. It is easy to intuitively understand what it is in the 1D case. Suppose you have a probability distribution for some random variable $x$ parametrized by $\theta$: $P(x|\theta)$. If you change $\theta$ by an infinitesimal amount, how will the probability distribution itself change? Will it be vastly different or almost the same? We can quantify that change by averaging the square of the relative changes of the probabilities of all the points: $$\mathcal{I}=E\left[\left(\frac{dP(x)}{d\theta}\frac{1}{P(x)}\right)^2\right]_x$$

Another way to look at it is as quantifying the “resolution” with which we can detect the parameter $\theta$: when the FI is high, we can distinguish between 2 parameters with close values more easily than for low FI, which corresponds to a higher resolution in parameter space. But why can I make this statement? After all the FI quantifies the difference between the probability distributions, not the parameters which specify them. The reason the “resolution” picture makes sense is thanks to the Cramér-Rao bound: $$var(\hat{\theta})=\frac{1}{\mathcal{I}}$$

To understand the bound we must make the following definitions: a hat over a parameter denotes the estimator of that parameter (an estimator is just a function that takes a set of realizations (or “measurements”) of the random variable we are looking at and spits out a number, which we hope will be close to the true value of the estimator) and the variance of an estimator is just the variance of the result of applying the estimator to many independent sets of measurements. With that in mind, the statement is that the variance of an estimator of our parameter is bounded by the inverse of the FI (note: we assume the estimator to be “unbiased“). Because FI for a set of N independent samples scales as N, for a very large sample size the variance of the estimator tends to zero, meaning we always get the same, “correct value”, as we’d intuitively expect.

Now bear in mind this is a very mathematical construct, hence all the quotation marks. The intuition I described above is only valid when the model is simple enough that we can afford to use the concept of unbiased estimators (which imply discarding all prior information we might have) and can assume that the underlying parameter is indeed a single number. Im many interesting cases we cannot assume this, but rather that it is a random variable itself (either due to nature, or due to intermediate processes we neglect to add to our model). I digress. The FI stands on its own as a useful tool in a myriad of applications, often related to quantifying the “resolution” of a system in the sense I tried to convey above.

In the paper we are looking at, the authors consider dynamic systems at a microscopic level with many degrees of freedom. In these systems you can consider the attributes of each particle a parameter, and code the many possible state of the system at some point in time by a probability distribution, which of course depends on this very large set of parameters. In this case the FI for any parameter is likely very small, as if you were to make a small change in one of these parameters the system would evolve in a very similar way leaving the probability distributions relatively unchanged, in the FI sense. The interesting step is now to find the eigenvalues of the FI matrix. This projects the parameters onto a new space where the directions correspond to some natural observables of the system.

Now, if we coarsen the system by allowing a long time to pass (i.e. a diffusion process) or by looking at it from a macroscopic scale (i.e. coarse graining an ising model) it turns out a few of these directions have a very large weight and the rest have comparatively low weight. The authors argue that these directions, when cast as observables, correspond to the macroscopic parameters of the system. Going back to the picture of the FI as resolution, these few observables will be the ones which we will be able to easily distinguish, while all the others will get lost in the noise. This is an appealing statement because it agrees with what we already know from statistical physics: we can accurately model systems at the macroscopic scale even if we have no hope to know what is going on at the microscopic level. Now we can see this idea emerge naturally from probability theory.

Another point they make is that this procedure works for both diffusive type processes, where we attribute this scale separation due to the fact that fluctuations are only relevant at the micro scale but not at the macro; and for processes with phase transitions where fluctuations are relevant at all scales at the critical point (cf. renormalization group). Under this framework there is a single explanation for why universal behavior is so prevalent in physics which I think is pretty cool.

Modeling excitable media with cellular automata

While researching for my seminar I came across a class of cellular automata which models spiral waves in excitable media. Because these models are so simple I had some fun implementing them in processing. Processing is great because you can use the javascript version to embed the visualization in a webpage directly and everyone can play with it. Here are some of the models I played around with:

  • Spiral, a model developed by Gerhard and Schuster to simulate a chemical reaction.
  • CCA, a simple cyclic cellular automaton.
  • StochCCA, where I added stochasticity to the previous model. This is useful to assign more weight to the 4 neighbors of the Von Neumann neighborhood than to the remaining 4 which complete the Moore neighborhood. This appears to make space “more isotropic” and makes the waves actually circular.

Rays for android

Another platform I had a go at porting the rays app for was android. This was probably the least enjoyable port I’ve done, not because of the Java language but because there is essentially no decent library to do the kinds of 2D effects the app requires.

I researched a bit on the topic and there were three ways to go about doing this:

  • Use android’s own openGL implementation.
  • Use the NDK to write native code which draws directly to the framebuffer.
  • Use libGDX.

At the time (well over a year ago, before the fancy new android developers redesign) documentation for the first two options was absolutely incomprehensible. Even now I find using the openGL implementation in android quite confusing, although it might be my fault as I never really learned low level openGL. I now feel tempted to use the much improved NDK and do some low level pixel manipulation, but other projects take priority…

I ended up using libGDX. It was also poorly documented (some random wiki pages + code comments) but the API made a lot more sense to me. This is because it was in the vein of some 2D libraries for the NDS such as the venerable PAlib and uLibrary which I used extensively previously (in hindsight, it appears they were not the best libraries in terms of code quality and the simple libNDS would have been a much better choice). The fact that you could instantly play the game on your computer natively without needing an emulator also helped. Though I would recommend just plugging in your android device and debugging directly there, Android’s debugger is truly excellent.

The library was chosen, the code was banged out, et voila! It lives! No antialiasing, no funky glow effects but it lives!

Rays app running on Galaxy S2
Rays app running on Galaxy S2

One idea I had to make the app more interesting would be to synthesize some musical instrument, something like a violin, where each instrument would be represented by a ray, with the intensity of the sound proportional to its speed. I could not find any satisfactory samples though, since they all come as single notes and I’d rather have some continuous sound produced. It might be more interesting to use some sort of synthesis algorithm to produce an interesting sound and modulate its amplitude, frequency or harmonics with the rays. How to go about it remains a mystery  There aren’t too many synthesis libraries available for android that do what I want. Libpd seems to be a decent option, now I just need to find the time to finish this project.

Since the application has only been tested on a handful of phones, I never did toss it on the play store. Get it here if you want to play with it. You’ll have to allow installation of non market apps.

Simulating networks of nonlinear stochastic systems

arXiv:1209.3700

In this paper we attempt to find a computationally efficient way to numerically simulate networks with nonlinear stochastic dynamics. With this I mean a continuous dynamical model where the differential equation for each variable depends nonlinearly on some or all variables of the system and has additive noise. If $x$ is a vector with all variables and $\eta$ is a random vector of the same size as $x$ with some unspecified distribution, the dynamics can be compactly described as $$\frac{d x}{dt}=f(x,t)+\eta$$

The challenge lies in the nonlinearity combined with stochasticity. Were only one of them to be present, the problem would be simple. A deterministic nonlinear problem can be straightforwardly be integrated with an ODE package, while a linear stochastic system can be reduced to a system of ODEs for the moments of the probability distribution function (PDF). A full solution would require a Monte Carlo algorithm to simulate a sufficient number of paths to allow us to estimate the PDF of $x$ at each time point. For networks with many nodes we are haunted by the curse of dimensionality, as the volume needed to be sampled increases exponentially and so do the number of simulated paths required to get a good approximation of the distribution at later time points. In systems where there is a well defined mode around which most of the probability mass is concentrated we should be able to derive an analytic approximation which is more tractable. This is exactly what we try to do in the paper. Continue reading “Simulating networks of nonlinear stochastic systems”