Title: Dynamic Properties of Networks Author: Thomas M. Breuel Institution: UniKL
from pylab import *
(Examples of Networks)
In this chapter, we consider the dynamics of networks. Most of our examples will be from the dynamics of interacting neurons. However, the same phenomena occur in many other kinds of networks:
(Examples of Variables)
(Phenomena)
Common phenomena are:
These occur even if individual network components are non-oscillatory and deterministic, and in systems with both continuous and discrete time.
(Practical Significance)
These phenomena have enormous practical significance. For example:
We will primarily focus on networks of neurons today.
Oscillations are observed widely throughout the nervous system. We have previously seen oscillatory behavior in individual neurons. It turns out that oscillatory behavior can also arise in networks of neurons, even spontaneously.
This kind of oscillatory behavior is common in many other systems, including economic systems, computer networks, etc. It can be induced even by weak or sparse coupling between individual system components.
Let us start with a simple network of weakly coupled integrate-and-fire neurons. As you can see, although the neurons start out firing synchronously and firing causes positive feedback, the synchronicity is actually lost over time.
# weakly coupled integrate-and-fire neurons
N = 1000
state = 0.001*rand(N)
coupling = 0.001*rand(N,N)
activity = []
for trial in range(5000):
state += 0.01*rand(N)
output = (state>1.0)
state[output] = 0.0
state += dot(coupling,output)
activity.append(sum(output))
print amin(activity),amax(activity)
0 653
# weakly coupled integrate-and-fire neurons
figsize(12,4)
plot(activity)
[<matplotlib.lines.Line2D at 0x43c0ed0>]
We can achieve synchronicity by introducing a short refractory period after a neuron fires. This leads to global synchronization even if neurons start out in completely unsynchronized states.
N = 1000
state = rand(N)
tref = 5
refr = zeros(N)
coupling = 0.001*rand(N,N)
activity = []
for trial in range(5000):
state += 0.01*rand(N)
refr = maximum(0,refr-1)
output = (state>1.0)
refr[output] = tref
state[refr>0] = 0.0
state += dot(coupling,output)
activity.append(sum(output))
print amin(activity),amax(activity)
0 634
# weakly coupled integrate-and-fire neurons with refractory period
plot(activity)
[<matplotlib.lines.Line2D at 0x4509550>]
We can look at the shape of an individual peak. We see an exponential growth of neurons firing, as each firing triggers even more neurons to fire. Then, however, all neurons cease simultaneously and their internal state variables start from zero again.
plot(activity[4800:5000])
[<matplotlib.lines.Line2D at 0x49ee610>]
The refractory period allows neurons to be reset to zero after a peak has happened.
The refractory period needs to be about as long as the peak itself is wide.
If the refractory period is not present, neurons active during the peak will partially activate subsequent neurons, leaving the population with a mix of initial conditions.
If the refractory period is present, neurons active during the peak will only become sensitive again after most other neurons have stopped firing, so the active population will start with a similar phase value or internal state, causing synchronization.
(Remember)
Synchronization can occur as a result of weak interactions in a network and refractory periods.
Synchronization can also occur by weak coupling of non-linear oscillators.
In fact, with nonlinear oscillators, we can see a stronger phenomenon, namely entrainment, where oscillators of different frequencies adjust not just their phase but also their frequency.
(History)
Synchronization and entrainment were first observed with pendulum clocks.
They tend to have slightly different frequencies and different phases, but they would synchronize when placed on the same table.
Let's start with the van der Pol oscillator.
Remember that the van der Pol oscillator is a parameter choice for the FitzHugh-Nagumo model, and hence the phenomena described here also apply to interacting neurons.
The differential equations are:
$$ \dot{x} = y $$$$ \dot{y} = \mu (1-x^2) y - x $$Here, $\mu$ is a parameter that determines the form of the oscillations. For $\mu=0$, we get regular harmonic oscillations.
The effect off the non-linear damping is to add energy to the system if there is too little energy in the system, and to remove energy from the system if there is too much present.
# van der Pol oscillator
x,y = 2,-4 # 0.1,0.2
mu = 1.0
s = 0.01
l = []
for i in range(10000):
dx = y
dy = mu * (1-x**2) * y - x
x += dx*s
y += dy*s
l.append((x,y))
l = array(l)
# output from the van der Pol oscillator
plot(l[:,0])
[<matplotlib.lines.Line2D at 0x4d7ee50>]
As usual, we can plot this in phase space.
The first nullcline is trivial:
$$ \dot{x} = y = 0 $$The second nullcline is quite interesting:
$$ \dot{y} = \mu (1-x^2) y - x = 0 $$or
$$ y = \frac{x}{\mu (1-x^2)} $$# van der Pol oscillator in phase space
plot(l[:,0],l[:,1])
xs = linspace(-3,3,117)
plot(xs,0*xs)
plot(xs,clip(xs/(1-xs**2)/mu,-5,5))
[<matplotlib.lines.Line2D at 0x4a027d0>]
Let's look at two van der Pol oscillators running independently, with different starting points and different $\mu$ parameters (frequencies). As you can see, if these run independently, there is no fixed relationship between the two oscillators.
# uncoupled van der Pol oscillators
N = 2
x,y = randn(N),randn(N)
mu = 1+0.5*arange(N)
s = 0.1
l = []
for i in range(20000):
dx = y
dy = mu * (1-x**2) * y - x
x += dx*s
y += dy*s
l.append(x.copy())
l = array(l)
# uncoupled van der Pol oscillators
plot(l[-1000:,0])
plot(l[-1000:,1])
[<matplotlib.lines.Line2D at 0x4e0e950>]
We can actually plot one oscillator against the other. The dense filling of space tells us again that there is no fixed relationship.
# uncoupled van der Pol oscillators (trajectory)
plot(l[-3000:,0],l[-3000:,1])
[<matplotlib.lines.Line2D at 0x5c7cb50>]
Now we introduce a small amount of mutual coupling between the two. As you can see, this causes both the frequency and the phase of the two oscillators to align.
# coupled van der Pol oscillators
N = 2
x,y = 1.0+arange(N),0.0+arange(N)
mu = 1+1.5*arange(N)
s = 0.1
l = []
for i in range(10000):
dx = y
dy = mu * (1-x**2) * y - x + 0.4*roll(dx,1)
x += dx*s
y += dy*s
l.append(x.copy())
l = array(l)
# output from coupled van der Pol oscillators
plot(l[-1000:,0])
plot(l[-1000:,1])
[<matplotlib.lines.Line2D at 0x5c84510>]
When coupling leads to agreement of the frequency of two oscillators (at some fixed phase), it is called entrainment of the two oscillators. It is a common phenomenon in dynamical systems.
You can actually try to synchronize your brainwaves to an external source (or another person's brainwaves) in a process also called entrainment.
# entrainment in phase space
plot(l[:,0],l[:,1],c=(0,0,1,0.3))
[<matplotlib.lines.Line2D at 0x58c8dd0>]
(Summary)
Weak coupling of nonlinear oscillators can lead to entrainment, the synchronization of both phase and frequency of the two oscillators.
Above, the two van der Pol oscillators were coupled weakly and symmetrically.
We can also use one oscillator to drive the other oscillator strongly.
In that case we get chaotic behavior.
# forced van der Pol oscillator
N = 2
x,y = 1.0+arange(N),0.0+arange(N)
mu = 1+1.5*arange(N)
s = 0.1
l = []
for i in range(10000):
dx = y
dy = mu * (1-x**2) * y - x
dy += 2.0 * roll(dx,1)*(arange(N)>0) # COUPLING
x += dx*s
y += dy*s
l.append(x.copy())
l = array(l)
# chaotic behavior of forced van der Pol oscillator
plot(l[-1000:,0])
plot(l[-1000:,1])
[<matplotlib.lines.Line2D at 0x4e2f310>]
# forced oscillator in phase space
plot(l[-3000:,0],l[-3000:,1])
[<matplotlib.lines.Line2D at 0x65f4cd0>]
The forced van der Pol oscillator was the first system where deterministic chaos was observed.
Van der Pol had built an electronic oscillator out of vacuum tubes and listening to the signal.
When he tried to drive that oscillator with another oscillator, at certain parameter choices, noise was heard, instead of interference or "beating" between the oscillators.
(Summary)
Strong forcing of a non-linear oscillator can lead to chaotic dynamics.
Chaotic behavior is similar to true randomness.
The Lorentz system is another differential equation with chaotic behavior.
It occurs in weather and climate modeling, as well as many other areas of science and economics.
It has three time-dependent variables and is given by:
$$ \dot{x} = \sigma(y-x) $$$$ \dot{y} = x(\rho-z)-y $$$$ \dot{z} = xy - \beta z $$# Lorenz system
x,y,z = 1.0,2.0,0.0
sigma,beta,rho = 10.0,8.0/3,28.0
s = 0.001
l = []
for i in range(100000):
dx = sigma*(y-x)
dy = x*(rho-z)-y
dz = x*y-beta*z
x = clip(x+dx*s,-1e3,1e3)
y = clip(y+dy*s,-1e3,1e3)
z = clip(z+dz*s,-1e3,1e3)
l.append((x,y,z))
l = array(l)
# trajectory of Lorenz system
plot(l[-10000:,0])
[<matplotlib.lines.Line2D at 0x66236d0>]
The usual way of plotting this is in three-dimensional phase space.
Note that these phase space trajectories never converge to a single limit cycle as they do for other oscillators.
Instead, they move around in a lower dimensional subspace, a strange attractor.
from mpl_toolkits.mplot3d.axes3d import Axes3D
ax = gcf().add_subplot(1, 1, 1, projection='3d')
figsize(8,8)
ax.plot(l[:,0],l[:,1],l[:,2],color=(0,0,1,0.8))
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x6823510>]
Let's simulate multiple paths with slightly different initial conditions.
This shows sensitivity to initial conditions or the butterfly effect.
N = 10
x,y,z = ones(N)+1e-6*arange(N),zeros(N),zeros(N)
sigma,beta,rho = 10.0,8.0/3,28.0
s = 0.001
l = []
for i in range(100000):
dx = sigma*(y-x)
dy = x*(rho-z)-y
dz = x*y-beta*z
x = clip(x+dx*s,-1e3,1e3)
y = clip(y+dy*s,-1e3,1e3)
z = clip(z+dz*s,-1e3,1e3)
l.append(x.copy())
l = array(l)
With this simulation, initially, the trajectories are close together.
figsize(12,4)
plot(l[:30000,0]); plot(l[:30000,1]); plot(l[:30000,-1])
[<matplotlib.lines.Line2D at 0x84d03d0>]
But they soon start to diverge. Note that there is a point near the beginning where one path goes into a completely different direction from the others.
figsize(12,4)
t = 30000
plot(l[t:t+30000,0]); plot(l[t:t+30000,1]); plot(l[t:t+30000,-1])
[<matplotlib.lines.Line2D at 0x84fddd0>]
Eventually, the trajectories have all largely diverged.
figsize(12,4)
plot(l[-30000:,0]); plot(l[-30000:,1]); plot(l[-30000:,-1])
[<matplotlib.lines.Line2D at 0x7d5b390>]
(Summary)
The Lorenz equations describe a chaotic dynamical system.
These equations exhibit sensitive dependence to initial conditions, an important aspect of chaotic behavior.
(Note)
Not all sensitive dependence to initial conditions results in chaos.
For example, the equation $\dot{x} = r\;x$ also has sensitive dependence to initial conditions (leading to exponential growth), but the behavior is not chaotic.
The van der Pol oscillator and Lorentz equations both may seem a little obscure and complicated.
Even simpler systems can exhibit chaos. Of those, the Lotka-Volterra equations are the most important, since they describe the behavior of many kinds of systems in economics, social networks, e-commerce, etc.
# Lotka-Volterra Equations
x = 0.1
y = 0.1
s = 0.01
l = []
for i in range(10000):
dx = x * (1-y)
dy = -0.3 * y * (1-x)
x += dx * s
y += dy * s
l.append((x,y))
# Lotka-Volterra Output
l = array(l)
plot(l)
[<matplotlib.lines.Line2D at 0x8d6cb50>, <matplotlib.lines.Line2D at 0x8d6cdd0>]
We can extend the Lotka-Volterra equations to multiple species and make them competitive.
That is, the change in population $x$ is given by:
$$ \dot{x_i} = r \; x_i \; (1 - x_i - \sum \beta_j x_j) $$That is, in addition to the regular growth limiting factor $1-x_i$, growth is also limited by competition from other species.
# competitive Lotka-Volterra equations
N = 4
x = rand(N)
r = array([1,.72,1.53,1.27])
M = array([[1,1.09,1.52,0],[0,1,0.44,1.36],[2.33,0,1,0.47],[1.21,0.51,0.35,1]])
print M
s = 0.1
l = []
for i in range(50000):
dx = r * x * (1-dot(M,x))
x += dx * s
l.append(x.copy())
[[ 1. 1.09 1.52 0. ] [ 0. 1. 0.44 1.36] [ 2.33 0. 1. 0.47] [ 1.21 0.51 0.35 1. ]]
# competitive Lotka-Volterra equations behavior
l = array(l)
_=plot(l[-10000:])
We can also plot this in a higher dimensional phase space, and we see a strange attractor again. (We're starting the plot at $t=10000$ after the transients have settled.)
from mpl_toolkits.mplot3d.axes3d import Axes3D
ax = gcf().add_subplot(1, 1, 1, projection='3d')
figsize(8,8)
ax.plot(l[10000:,2],l[10000:,1],l[10000:,0],color=(0,0,1,0.8))
[<mpl_toolkits.mplot3d.art3d.Line3D at 0xc424910>]
(Summary)
Simple competitive growth (as occurs in economics, markets, ecology, and many other fields) can lead to chaotic behavior.
All the above models are differential equation models, or smooth, continuous dynamical systems.
(We simulate them as difference equations but choose the step size small enough to get continuous trajectories).
However, even simple difference equations can exhibit chaotic behavior.
These are some of the simplest dynamical models exhibiting chaos.
(Ricker's Map)
Ricker's map simulates bounded growth:
$$ x_{t+1} = x_t \; e^{r(1-x_t/k)} $$The growth rate is given by $ e^{r(1-x/k)} - 1$.
Once the population goes above $k$, growth rates become negative and the population declines
# Ricker's map
xs = []
x = 0.1
k = 1.0
r = 3.0
for i in range(100):
xs.append(x)
x = x * exp(r*(1-x/k))
plot(xs)
[<matplotlib.lines.Line2D at 0x8e70050>]
(Logistic Map)
An even simpler example is the logistic map, a discrete analog of the Lotka-Volterra equations and purely a polynomial:
$$ x_{t+1} = r \; x_t \; (1-x_t) $$# Logistic map
xs = []
x = 0.1
r = 3.9
for i in range(100):
x = r * x * (1-x)
xs.append(x)
plot(xs)
[<matplotlib.lines.Line2D at 0x9539710>]
(Summary)
Chaotic behavior can occur even in simple discrete time systems and with polynomial difference equations like the logistic map.
There is another system where chaotic behavior can be observed, namely cellular automata.
Chaos appears even in the simplest kind of cellular automata.
# Generic Rule N implementation
def rule_n(a,n):
m = array([(n&(1<<i)!=0) for i in range(8)],'i')
s = array(4*roll(a,1)+2*a+roll(a,-1),'i')
return m[s]
# run a 1D cellular automaton through several iterations
def run(n,initial=None):
if initial is not None:
a = initial
else:
a = zeros(500)
a[350] = 1
result = []
for i in range(400):
result.append(a.copy())
a = rule_n(a,n)
return array(result)
# Rule 110 is at the edge of chaos
gray(); figsize(12,12)
imshow(run(110))
<matplotlib.image.AxesImage at 0xa20c110>
(Rule 110 and the Border of Chaos)
# Rule 30 is chaotic
gray(); figsize(12,12)
imshow(run(30))
<matplotlib.image.AxesImage at 0xa232ed0>
(Rule 30 and Chaos)
(Classification of Cellular Automata)
(Classification of Cellular Automata)
(Kinds of Behavior)
(Main Message)
Even random networks of simple components with simple interactions can exhibit much more complex behaviors than the individual components.
These kind of behaviors are called emergent behaviors.
These complex phenomena occur in many areas of computer science, social systems, and economics.
(Questions)
Think about the kinds of software systems you create.
Where do you make assumptions about average cases or statistical independence?
Where is there a potential for weak coupling?
Where is there a potential for chaos?