In [1]:
%matplotlib widget
from tvb.simulator.lab import *
from tvb.simulator.plot.phase_plane_interactive import PhasePlaneInteractive

Exploring A Model

This tutorial discusses using the phase_plane_interactive plotting tool to explore the dynamics of a Model object and, at the same time, set its parameters.

This works best for the simpler, 2d, Models, as their trajectories and nullclines lie within a single plane and are thus easily visualised. And so, for the demo here we'll stick with a Model of that type. However, although it requires a little more effort, it can still be used to get a basic handle on the dynamics of higher dimensional models.

It is also important to note that this is only for the local dynamic model, that is, it only represents the dynamic behaviour of a disconnected node.

Create An Interactive Phase-Plane

PhasePlaneInteractive produces an interactive window containing plots of a Model's phase-plane, sample trajectories, and sliders and buttons for adjusting various parameters.

The important thing to note here is that as we drag around the sliders for the Model's parameters we are actually modifying the parameters of the Model we passed in, so when we close the figure at the end whatever values the Model's parameters are set to will be the values in the Model. Also, changing the range of the phase-plane plot (that is, the extent of the x and y axis) changes the "state_variable_range" attribute of a Model. This attribute is used when constructing random initial conditions for a simulation, so setting the axis to be relatively tightly bound around a fixed point for example will produce initial conditions that better approximate the range of the Model's state variables for the given parameters.

We'll begin by creating a Model object and taking a quick look at its parameters:

In [2]:
oscillator = models.Generic2dOscillator()

oscillator
Out[2]:

Generic2dOscillator

value
I [min, median, max]
[0, 0, 0]
I dtype
float64
I shape
(1,)
Type
Generic2dOscillator
a [min, median, max]
[-2, -2, -2]
a dtype
float64
a shape
(1,)
alpha [min, median, max]
[1, 1, 1]
alpha dtype
float64
alpha shape
(1,)
b [min, median, max]
[-10, -10, -10]
b dtype
float64
b shape
(1,)
beta [min, median, max]
[1, 1, 1]
beta dtype
float64
beta shape
(1,)
c [min, median, max]
[0, 0, 0]
c dtype
float64
c shape
(1,)
d [min, median, max]
[0.02, 0.02, 0.02]
d dtype
float64
d shape
(1,)
e [min, median, max]
[3, 3, 3]
e dtype
float64
e shape
(1,)
f [min, median, max]
[1, 1, 1]
f dtype
float64
f shape
(1,)
g [min, median, max]
[0, 0, 0]
g dtype
float64
g shape
(1,)
gamma [min, median, max]
[1, 1, 1]
gamma dtype
float64
gamma shape
(1,)
gid
UUID('1cc8c9ac-96cf-40fa-85e5-ea046d932969')
state_variable_range
{'V': array([-2.,  4.]), 'W': array([-6.,  6.])}
tau [min, median, max]
[1, 1, 1]
tau dtype
float64
tau shape
(1,)
title
Generic2dOscillator gid: 1cc8c9ac-96cf-40fa-85e5-ea046d932969
variables_of_interest
('V',)

We'll now create and launch the interactive plot.

NOTE: Changing the Model's parameters or the axis settings causes a redraw of the entire phase-plane, clearing trajectories and their corresponding time-series.

In [3]:
ppi_fig = PhasePlaneInteractive(model=oscillator)
ppi_fig.show()

In the main central panel of the window you can see the phase-plane for the model, including arrows representing the vector field and coloured lines representing any nullclines present in this plane. Clicking on the phase-plane will launch a sample trajectory originating from where you clicked. Below the phase-plane is a panel which will show time-series for all state variables for any sample trajectories you initiate. All around the edges are sliders for adjusting Model parameters and adjusting what is displayed. The red vertical lines in sliders indicate the initial values.

After we've adjusted parameters to our satisfaction we can close the window and take another quick look at the parameters of our Model

In [4]:
oscillator
Out[4]:

Generic2dOscillator

value
I [min, median, max]
[0, 0, 0]
I dtype
float64
I shape
(1,)
Type
Generic2dOscillator
a [min, median, max]
[-2, -2, -2]
a dtype
float64
a shape
(1,)
alpha [min, median, max]
[1, 1, 1]
alpha dtype
float64
alpha shape
(1,)
b [min, median, max]
[-10, -10, -10]
b dtype
float64
b shape
(1,)
beta [min, median, max]
[1, 1, 1]
beta dtype
float64
beta shape
(1,)
c [min, median, max]
[0, 0, 0]
c dtype
float64
c shape
(1,)
d [min, median, max]
[0.02, 0.02, 0.02]
d dtype
float64
d shape
(1,)
e [min, median, max]
[3, 3, 3]
e dtype
float64
e shape
(1,)
f [min, median, max]
[1, 1, 1]
f dtype
float64
f shape
(1,)
g [min, median, max]
[0, 0, 0]
g dtype
float64
g shape
(1,)
gamma [min, median, max]
[1, 1, 1]
gamma dtype
float64
gamma shape
(1,)
gid
UUID('1cc8c9ac-96cf-40fa-85e5-ea046d932969')
state_variable_range
{'V': array([-2.,  4.]), 'W': array([-6.,  6.])}
tau [min, median, max]
[1, 1, 1]
tau dtype
float64
tau shape
(1,)
title
Generic2dOscillator gid: 1cc8c9ac-96cf-40fa-85e5-ea046d932969
variables_of_interest
('V',)

As you can see in the line above, the Model's parameters, for example "a", "tau", and the state_variable_ranges are modified from their initial values.

Specifying Stochastic Integration

It is possible to explicitly specify the integration scheme used to plot sample trajectories. This can be useful when deciding what amplitude to give your noise when specifying a stochastic integration scheme.

We'll take a look at this using HeunStochastic, we'll also pass in the same Model object we modified above. In this way PhasePlaneInteractive initialises with the parameters we'd set for the Model, so that here we can focus on the effect of the noise amplitude relative to the intrinsic dynamics.

Unlike changes to Model parameters and the axes, changing to the noise amplitude doesn't cause a redraw of the existing trajectories, so, after creating a trajectory you can alter the noise strength and click on the same starting location to see the effect of a different noise amplitude on the same trajectory.

Starting by setting the noise to 0.0, to get a deterministic trajectory, and then adding a small amount of noise can help give a useful intuition for the effects of noise on a simulation.

Also, as the random sequence used each time you launch a trajectory is distinct, clicking on the same point multiple times will give you an idea of the range of trajectory a given type of noise can produce.

Alternatively, clicking reset random stream enables you to see the effect of the same random sequence with different amplitudes or on a trajectory initiating from a different location.

Note: The default standard deviation for the noise is 1.0, which is too large (in the sense that noise will dominate the intrinsic dynamics) relative to the range of our Model's state variables.

In [5]:
import numpy as np

nsigma = np.array([0.01,])
hiss = noise.Additive(nsig=nsigma)
heunstochint = integrators.HeunStochastic(dt=2**-5, noise=hiss)

# examine integrator's attributes
heunstochint
Out[5]:

HeunStochastic

value
Type
HeunStochastic
bounded_state_variable_indices
None
clamped_state_variable_indices
None
clamped_state_variable_values
None
dt
0.03125
gid
UUID('1e935fd3-7cf0-403f-9e0e-4ef424d97f48')
noise
Additive gid: 6b62573b-a5b2-4658-b3bb-3a6941e35c63
state_variable_boundaries
None
title
HeunStochastic gid: 1e935fd3-7cf0-403f-9e0e-4ef424d97f48
In [6]:
heunstochint.noise
Out[6]:

Additive

value
Type
Additive
gid
UUID('6b62573b-a5b2-4658-b3bb-3a6941e35c63')
noise_seed
42
nsig [min, median, max]
[0.01, 0.01, 0.01]
nsig dtype
float64
nsig shape
(1,)
ntau
0.0
random_stream
RandomState(MT19937) at 0x7F7C90D6E740
title
Additive gid: 6b62573b-a5b2-4658-b3bb-3a6941e35c63

Relaunch the phase plane tool, but with the stochastic integrator

In [8]:
ppi_fig = PhasePlaneInteractive(model=oscillator, integrator=heunstochint)
ppi_fig.show()
In [9]:
heunstochint.noise
Out[9]:

Additive

value
Type
Additive
gid
UUID('6b62573b-a5b2-4658-b3bb-3a6941e35c63')
noise_seed
42
nsig [min, median, max]
[0.00471795, 0.00471795, 0.00471795]
nsig dtype
float64
nsig shape
(1,)
ntau
0.0
random_stream
RandomState(MT19937) at 0x7F7C90D6E740
title
Additive gid: 6b62573b-a5b2-4658-b3bb-3a6941e35c63

Simulate and Visualize

Finally, we can use the objects created above in a simulation:

In [12]:
conn = connectivity.Connectivity.from_file()
period = 2**-1
sim = simulator.Simulator(
    model = oscillator, 
    connectivity = conn,
    coupling = coupling.Linear(a=np.array([0.0152])), 
    integrator = heunstochint, 
    monitors = (monitors.TemporalAverage(period=period),),
    simulation_length=1e3,
).configure()

# run
(tavg_time, tavg_data), = sim.run()
In [15]:
#Create a TVB TimeSeries object
import tvb.datatypes.time_series
tsr = tvb.datatypes.time_series.TimeSeriesRegion()
tsr.data = tavg_data
tsr.sample_period = period /1000
tsr.sample_period_unit = 's'
tsr.connectivity = conn

from tvb.simulator.plot.timeseries_interactive import TimeSeriesInteractivePlotter
tsi = TimeSeriesInteractivePlotter(time_series=tsr)
tsi.configure()
tsi.show()

As long as the interactive phase plane figure isn't closed, the parameters can be tuned and the simulation rerun to iterate between the local dynamics of a node and integration scheme and the global dynamics of the entire network.

In [ ]: