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
)
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
)
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
)
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
)
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
)
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
)
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
)