Continuing a region stochastic BOLD simulation

In many cases, it is useful to perform an initial simulation to allow the transient dynamics to dissipate, and then continue the simulation using the steady state obtained from the initial simulation as the initial conditions. The TVB web interface allows for this by 'branching' a loaded simulation. From a Python script, we need to handle the details ourselves. This demo shows how to do that for a region-level, stochastic simulation with the BOLD monitor.

In [12]:
%pylab inline
from tvb.simulator.lab import *
LOG = get_logger('demo')
import cPickle
Populating the interactive namespace from numpy and matplotlib

Initial transient

First, we build and configure the initial simulator,

In [13]:
sim = simulator.Simulator(
    model=models.Generic2dOscillator(), 
    connectivity=connectivity.Connectivity(load_default=True),                      
    coupling=coupling.Linear(a=0.0126),
    simulation_length=1e4,
    integrator=integrators.HeunStochastic(
        dt=2 ** -4,
        noise=noise.Additive(nsig=array([0.001]))),
    monitors=(
        monitors.TemporalAverage(period=1.0),
        monitors.Bold(period=500.0),
        monitors.ProgressLogger(period=1e3)
    )
).configure()
WARNING  File 'hemispheres' not found in ZIP.

let it run,

In [14]:
(tavg_time, tavg_data), (bold_time, bold_data), _ = sim.run()
   INFO  step 1 time 0.0001 s
   INFO  step 16001 time 1.0001 s
   INFO  step 32001 time 2.0001 s
   INFO  step 48001 time 3.0001 s
   INFO  step 64001 time 4.0001 s
   INFO  step 80001 time 5.0001 s
   INFO  step 96001 time 6.0001 s
   INFO  step 112001 time 7.0001 s
   INFO  step 128001 time 8.0001 s
   INFO  step 144001 time 9.0001 s

and plot the transient:

In [15]:
figure(figsize=(10, 5))
plot(bold_time * 1e-3, bold_data[:, 0, :, 0], 'k', alpha=0.3)
title('Initial transient in BOLD')
xlabel('Time (s)')
grid(True);

To save the state of the simulator, we need to put a few arrays into files:

In [16]:
sim_state_fname = 'sim_state.pickle'

with open(sim_state_fname, 'wb') as file_descr:
    cPickle.dump({
        'history': sim.history.buffer,
        'current_step': sim.current_step,
        'current_state': sim.current_state,
        'bold_inner': sim.monitors[1]._interim_stock,
        'bold': sim.monitors[1]._stock,
        'rng': sim.integrator.noise.random_stream.get_state()
    }, file_descr)

The simulation state is now saved in the current folder:

In [17]:
!ls -lh sim_state.pickle
'ls' is not recognized as an internal or external command,
operable program or batch file.

We are free to dispose of the simulator:

In [18]:
del sim

Continuing the simulation

Now we want to continue the simulator above with the saved state. In different scenarios, this might be in a different script, on a different day, so we need to build the simulator object again:

In [19]:
sim = simulator.Simulator(
    model=models.Generic2dOscillator(), 
    connectivity=connectivity.Connectivity(load_default=True),                      
    coupling=coupling.Linear(a=0.0126),
    simulation_length=1e4,
    integrator=integrators.HeunStochastic(
        dt=2 ** -4,
        noise=noise.Additive(nsig=array([0.001]))),
    monitors=(
        monitors.TemporalAverage(period=1.0),
        monitors.Bold(period=500.0),
        monitors.ProgressLogger(period=1e3)
    )
).configure()
WARNING  File 'hemispheres' not found in ZIP.

load its state,

In [20]:
with open(sim_state_fname, 'rb') as file_descr:
    state = cPickle.load(file_descr)

sim.history.buffer = state['history']
sim.current_step = state['current_step']
sim.current_state = state['current_state']
sim.monitors[1]._interim_stock = state['bold_inner']
sim.monitors[1]._stock = state['bold']
sim.integrator.noise.random_stream.set_state(state['rng'])

and run it again; note that the step and time start off where the old simulator stopped:

In [21]:
(tavg_time, tavg_data), (bold_time, bold_data), _ = sim.run()
   INFO  step 160000 time 10.0000 s
   INFO  step 176000 time 11.0000 s
   INFO  step 192000 time 12.0000 s
   INFO  step 208000 time 13.0000 s
   INFO  step 224000 time 14.0000 s
   INFO  step 240000 time 15.0000 s
   INFO  step 256000 time 16.0000 s
   INFO  step 272000 time 17.0000 s
   INFO  step 288000 time 18.0000 s
   INFO  step 304000 time 19.0000 s
In [22]:
figure(figsize=(10, 5))
plot(bold_time * 1e-3, bold_data[:, 0, :, 0], 'k', alpha=0.3)
title('Continued BOLD')
xlabel('Time (s)')
grid(True);

The continued BOLD does not show the transient present in the initial simulation.