dynpy tutorial

Introduction

dynpy is a package for defining and running dynamical systems in Python. The goal is to support a wide-variety of dynamical systems, both continuous and discrete time as well as continuous and discrete state.

dynpy is organized into a hierarchy of classes, with each class representing a different kind of dynamical system. The base classes are defined in dynpy.dynsys - Base module for dynamical systems. Some definitions used for creating sample systems in this tutorial are defined in dynpy.sample_nets - Module with sample networks and dynamical systems.

Example: Boolean Networks

Boolean network are a type of discrete-state, discrete-time dynamical system. Each node updates itself as a Boolean function of other nodes (its ‘inputs’).

dynpy.bn - Boolean network module implement Boolean network dynamical systems. For example, we can use it to compute the space time diagram of the 11-node yeast cell-cycle network, as described in: Li et al, The yeast cell-cycle network is robustly designed, PNAS, 2004.

import matplotlib.pyplot as plt
import numpy as np
import dynpy

bn = dynpy.bn.BooleanNetwork(rules=dynpy.sample_nets.budding_yeast_bn)

initState = np.zeros(bn.num_vars, 'uint8')
initState[ [1,3,6] ] = 1
plt.spy(bn.get_trajectory(start_state=initState, max_time=15))
plt.xlabel('Node')
plt.ylabel('Time')

(Source code, png, hires.png, pdf)

_images/bn_trajectory.png

We can also get the network’s attractors, by doing:

>>> import dynpy
>>> bn = dynpy.bn.BooleanNetwork(rules=dynpy.sample_nets.budding_yeast_bn)
>>> atts, attbasins = bn.get_attractor_basins(sort=True)
>>> print(list(map(len, attbasins)))
[1764, 151, 109, 9, 7, 7, 1]

Or print them out using:

>>> import dynpy
>>> bn = dynpy.bn.BooleanNetwork(rules=dynpy.sample_nets.budding_yeast_bn)
>>> bn.print_attractor_basins()
* BASIN 1 : 1764 States
ATTRACTORS:
   Cln3    MBF    SBF Cln1,2   Sic1   Swi5  Cdc20 Clb5,6   Cdh1 Clb1,2   Mcm1
      0      0      0      0      1      0      0      0      1      0      0
--------------------------------------------------------------------------------
* BASIN 2 : 151 States
ATTRACTORS:
   Cln3    MBF    SBF Cln1,2   Sic1   Swi5  Cdc20 Clb5,6   Cdh1 Clb1,2   Mcm1
      0      0      1      1      0      0      0      0      0      0      0
--------------------------------------------------------------------------------
* BASIN 3 : 109 States
ATTRACTORS:
   Cln3    MBF    SBF Cln1,2   Sic1   Swi5  Cdc20 Clb5,6   Cdh1 Clb1,2   Mcm1
      0      1      0      0      1      0      0      0      1      0      0
--------------------------------------------------------------------------------
* BASIN 4 : 9 States
ATTRACTORS:
   Cln3    MBF    SBF Cln1,2   Sic1   Swi5  Cdc20 Clb5,6   Cdh1 Clb1,2   Mcm1
      0      0      0      0      1      0      0      0      0      0      0
--------------------------------------------------------------------------------
* BASIN 5 : 7 States
ATTRACTORS:
   Cln3    MBF    SBF Cln1,2   Sic1   Swi5  Cdc20 Clb5,6   Cdh1 Clb1,2   Mcm1
      0      0      0      0      0      0      0      0      0      0      0
--------------------------------------------------------------------------------
* BASIN 6 : 7 States
ATTRACTORS:
   Cln3    MBF    SBF Cln1,2   Sic1   Swi5  Cdc20 Clb5,6   Cdh1 Clb1,2   Mcm1
      0      1      0      0      1      0      0      0      0      0      0
--------------------------------------------------------------------------------
* BASIN 7 : 1 States
ATTRACTORS:
   Cln3    MBF    SBF Cln1,2   Sic1   Swi5  Cdc20 Clb5,6   Cdh1 Clb1,2   Mcm1
      0      0      0      0      0      0      0      0      1      0      0
--------------------------------------------------------------------------------

Cellular Automata

The cellular automata class dynpy.ca.CellularAutomaton is defined in dynpy.ca - Cellular automata module. It is a subclass of dynpy.bn.BooleanNetwork. It constructs a Boolean network with a lattice connectivity topology and a homogenous update function. For example:

import matplotlib.pyplot as plt
import numpy as np
import dynpy

ca = dynpy.ca.CellularAutomaton(num_vars=100, num_neighbors=1, rule=110)

initState = np.zeros(ca.num_vars, 'uint8')
initState[int(ca.num_vars/2)] = 1
plt.spy(ca.get_trajectory(start_state=initState, max_time=50))
plt.xlabel('Node')
plt.ylabel('Time')

(Source code, png, hires.png, pdf)

_images/ca_trajectory.png

Markov Chains

dynpy also implements Markov chains, or dynamical systems over distributions of states. See the documentation for dynpy.markov.MarkovChain for more details.

For example, here we use dynpy.graphdynamics - Module for implementing dynamical systems on graphs, which implements dynamics on graphs, to instantiate a dynamical system representing the distribution of a random walker on Zachary’s karate club network. Here, dynpy.graphdynamics.RandomWalker is a subclass of dynpy.markov.MarkovChain.

import matplotlib.pyplot as plt
import numpy as np
import dynpy

G = dynpy.sample_nets.karateclub_net
N = G.shape[0]
rw = dynpy.graphdynamics.RandomWalkerEnsemble(graph=G)

initState = np.zeros(N)
initState[ 5 ] = 1

trajectory = rw.get_trajectory(start_state=initState, max_time=30)
plt.imshow(trajectory, interpolation='none')
plt.xlabel('Node')
plt.ylabel('Time')

(Source code, png, hires.png, pdf)

_images/mc_trajectory.png

A Markov chain, like some other dynamical systems implemented by dynpy, can also be run in continuous time (in this context, it is sometimes called a ‘master equation’). This is specified by passing in the discrete_time=False option when constructing the underlying dynamical system. Using the previous example:

import matplotlib.pyplot as plt
import numpy as np
import dynpy

G = dynpy.sample_nets.karateclub_net
N = G.shape[0]
rw = dynpy.graphdynamics.RandomWalkerEnsemble(graph=G, discrete_time=False)

initState = np.zeros(N)
initState[ 5 ] = 1

trajectory = rw.get_trajectory(start_state=initState, max_time=30)
plt.imshow(trajectory, interpolation='none')
plt.xlabel('Node')
plt.ylabel('Time')

(Source code, png, hires.png, pdf)

_images/mc_cont_trajectory.png

It is also possible to get the equilibrium distribution by calling get_equilibrium_distribution(), which uses eigenspace decomposition:

import matplotlib.pyplot as plt
import numpy as np
import dynpy

kc = dynpy.sample_nets.karateclub_net
rw = dynpy.graphdynamics.RandomWalkerEnsemble(graph=kc, discrete_time=False)

eq_state = rw.get_equilibrium_distribution()
plt.imshow(np.atleast_2d(eq_state), interpolation='none')

(Source code, png, hires.png, pdf)

_images/equil_dist.png

In fact, it is possible to turn any deterministic dynamical system into a Markov chain by using the dynpy.markov.MarkovChain.from_deterministic_system() method. For example, to create a dynamical system over a distribution of states of the yeast-cell cycle Boolean network:

import matplotlib.pyplot as plt
import dynpy

bn = dynpy.bn.BooleanNetwork(rules=dynpy.sample_nets.budding_yeast_bn)
bnMC = dynpy.markov.MarkovChain.from_deterministic_system(bn)

# get distribution over states at various timepoints
t = bnMC.get_trajectory(start_state=bnMC.get_uniform_distribution(), max_time=20)

# project back from states onto activations of original nodes
bnProbs = t.dot(bn.get_ndx2state_mx())

# plot
plt.imshow(bnProbs, interpolation='none')
plt.xlabel('Node')
plt.ylabel('Time')

(Source code, png, hires.png, pdf)

_images/bn_dist.png

Stochastic Systems

Stochastic systems can also be implemented.

import matplotlib.pyplot as plt
import numpy as np
import dynpy

num_steps = 30
G = dynpy.sample_nets.karateclub_net
N = G.shape[0]
rw = dynpy.graphdynamics.RandomWalkerEnsemble(graph=G)
sampler = dynpy.markov.MarkovChainSampler(rw)

# Initialize with a single random walker on node id=5
trajectory = sampler.get_trajectory(start_state=5, max_time=num_steps)

plt.plot(np.arange(num_steps), trajectory, 'o')
plt.ylim([0, rw.num_states])
plt.xlabel('Time')
plt.ylabel('State')

(Source code, png, hires.png, pdf)

_images/random_walker.png