RuntimeWarning in simulate_region_Jansen_Rit

[Toygun Başaklar]

Hello,

I am working on Jansen Rit model to observe its response to visually evoked potentials. For this purpose, I have started with simulate_region_Jansen_Rit tutorial to practice using TVB. I encountered some RuntimeWarnings when I ran the code. Below I give the code and the corresponding warnings. How can I get rid of these warnings?

%pylab nbagg
from tvb.simulator.lab import *
jrm = models.JansenRit(mu=0., v0=6.)

phi_n_scaling = (jrm.a * jrm.A * (jrm.p_max-jrm.p_min) * 0.5 )**2 / 2.
sigma         = numpy.zeros(6)
sigma[3]      = phi_n_scaling

# the other aspects of the simulator are standard
sim = simulator.Simulator(
    model=jrm,
    connectivity=connectivity.Connectivity(load_default=True),
    coupling=coupling.SigmoidalJansenRit(a=10.0),
    integrator=integrators.HeunStochastic(dt=2 ** -4, noise=noise.Additive(nsig=sigma)),
    monitors=monitors.TemporalAverage(period=2 ** -2),
    simulation_length=1e3,
).configure()

# run it
(time, data), = sim.run()

# visualize time series
figure()
plot(time, data[:, 0, :, 0], 'k', alpha=0.1)
title("Temporal Average")

 
The output is as follows:

WARNING  File 'hemispheres' not found in ZIP.

/home/ee213/Desktop/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb/simulator/coupling.py:378: RuntimeWarning: overflow encountered in exp
  pre = self.cmax / (1.0 + numpy.exp(self.r * (self.midpoint - (x_j[:, 0] - x_j[:, 1]))))
/home/ee213/Desktop/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb/simulator/models/jansen_rit.py:310: RuntimeWarning: overflow encountered in _numba_dfun_jr
  self.A, self.b, self.B, self.J, self.mu

[Marmaduke Woodman]

Le 26 juil. 2016 à 09:35, Toygun Başaklar <tbas...@gmail.com> a écrit :

Hello,

I am working on Jansen Rit model to observe its response to visually evoked potentials. For this purpose, I have started with simulate_region_Jansen_Rit tutorial to practice using TVB. I encountered some RuntimeWarnings when I ran the code. Below I give the code and the corresponding warnings. How can I get rid of these warnings?

WARNING  File 'hemispheres' not found in ZIP.

This warning refers to an optional file, which isn’t present in the connectivity dataset you’re using. Again, it’s optional, so you can ignore this. 

/home/ee213/Desktop/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb/simulator/coupling.py:378: RuntimeWarning: overflow encountered in exp
  pre = self.cmax / (1.0 + numpy.exp(self.r * (self.midpoint - (x_j[:, 0] - x_j[:, 1]))))
/home/ee213/Desktop/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb/simulator/models/jansen_rit.py:310: RuntimeWarning: overflow encountered in _numba_dfun_jr
  self.A, self.b, self.B, self.J, self.mu

This exp term is part of a sigmoid of the form 1/(1 + exp(x)). Here, x is too large, exp(x) is, in floating point numbers, infinity. The result is still well defined, namely zero, so this warning can also be ignored.

To answer your actual question, such warnings are difficult to turn off, simply because they are often useful for debugging and not usually a problem to ignore. If you absolutely need to turn them off, please look up the warnings module in Python and similar topic in NumPy.

Cheers,

Marmaduke

[Toygun Başaklar]

Thank you for the quick reply. As far as I understand the second warning about _numba_dfun_jr is a follow up of the first warning, so they are practically the same. I still have some questions regarding the model:

1- So what you are saying is these warnings do not mean the simulation is wrong, am I right? I have implemented Jansen-Rit model in Matlab myself and I would like to compare its response with TVB's response to a flash stimulus, so I want to absolutely make sure this is working as intended.

2- Do I need to generate a stimulus as shown below, or can I add a transient exponential term to the "background activity noise" to account for stimulus as stated in the original paper?

eqn_t = equations.PulseTrain()
eqn_t.parameters["onset"] = 500.0 # ms
eqn_t.parameters["tau"]   = 5.0   # ms
eqn_t.parameters["T"]     = 500.  # 0.002kHz repetition frequency

stimulus = patterns.StimuliRegion(temporal = eqn_t,
                                  connectivity = conn,
                                  weight = stim_weights)

3- On a side note, Jansen's 1995 paper thoroughly discusses about a parameter K. I did not find any parameters related to this in JansenRit class. Does the connectivity matrix account for this K parameter?

[Marmaduke Woodman]

Le 26 juil. 2016 à 10:39, Toygun Başaklar <tbas...@gmail.com> a écrit :

1- So what you are saying is these warnings do not mean the simulation is wrong, am I right? I have implemented Jansen-Rit model in Matlab myself and I would like to compare its response with TVB's response to a flash stimulus, so I want to absolutely make sure this is working as intended. 

No, it doesn’t mean that the simulation is wrong.

2- Do I need to generate a stimulus as shown below, or can I add a transient exponential term to the "background activity noise" to account for stimulus as stated in the original paper?

Looks ok to me. You can always plot the trace of the stimulus to check it’s doing what you want

eqn_t = equations.PulseTrain()
eqn_t.parameters["onset"] = 500.0 # ms
eqn_t.parameters["tau"]   = 5.0   # ms
eqn_t.parameters["T"]     = 500.  # 0.002kHz repetition frequency

stimulus = patterns.StimuliRegion(temporal = eqn_t,
                                  connectivity = conn, 
                                  weight = stim_weights)

3- On a side note, Jansen's 1995 paper thoroughly discusses about a parameter K. I did not find any parameters related to this in JansenRit class. Does the connectivity matrix account for this K parameter?

The K parameter enters when JR consider two columns, so yes, in TVB, K is considered connectivity weights, not as parameters of the mass model itself. 

[Toygun Başaklar]

Okay thank you for clarifying the matter. I will continue writing on this topic if I have further inquiries about Jansen Rit Model.

--
You received this message because you are subscribed to a topic in the Google Groups "TVB Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/tvb-users/PWbmepE4zws/unsubscribe.
To unsubscribe from this group and all its topics, send an email to tvb-users+...@googlegroups.com.
To post to this group, send email to tvb-...@googlegroups.com.
Visit this group at https://groups.google.com/group/tvb-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/3C6E1D9B-4143-4337-B6F7-9C673B7ECF99%40gmail.com.

For more options, visit https://groups.google.com/d/optout.

[Toygun Başaklar]

I have a few more questions about the code and TVB in general. I'd be glad if you answered them.

1- Why does sigma have a size of 6? And how do we know 4th index corresponds to model input? Please do correct me if Im wrong.

2- Is there a way to visualize a noise or stimulus object? I couldn't find a number array to plot.


Thanks.

[Marmaduke Woodman]

Le 26 juil. 2016 à 14:11, Toygun Başaklar <tbas...@gmail.com> a écrit :

1- Why does sigma have a size of 6? And how do we know 4th index corresponds to model input? Please do correct me if Im wrong.

Sigma is the level of noise, and this model requires per state variable dispersion values. (so, it’s not the input, technically speaking)

2- Is there a way to visualize a noise or stimulus object? I couldn't find a number array to plot.

There’s a demo of visualizing the stimulus pattern, once configured. Please take a look and see if it helps.

Cheers,

Marmaduke

[Toygun Başaklar]

Sigma is the level of noise, and this model requires per state variable dispersion values. (so, it’s not the input, technically speaking)

But in terms of JR 1995, it is considered as the input to the system, isn't it? ("p(t)" in the below figure is sigma[3] in the code as far as I understood). I'd like to stimulate the model as described in the paper itself. In order to do that, I should add the desired stimulus (flash pulse in this case) to the noise but since the noise is not a data array (it is a noise object), I cannot add the stimulus. Can we somehow obtain access to this noise object's content (time,data) and manipulate it? 

 Figure from Jansen, Ben H., and Vincent G. Rit. "Electroencephalogram and Visual Evoked Potential Generation in a Mathematical Model of Coupled Cortical Columns." Biol. Cybern. Biological Cybernetics 73.4 (1995): 357-66. Web.

 
Another question that appeared to my mind, is applying a regional stimulus to this model in TVB equivalent to JR's way of applying stimulus (adding exponential onto white noise) in their paper?

There’s a demo of visualizing the stimulus pattern, once configured. Please take a look and see if it helps.

 

I took a look both at the tutorial and at the source code (tools.py). For stimulus it's perfectly explained but for noise I didn't find any relevant function to plot/visualize. The reason why I am trying to visualize the noise is that I want to be sure that I give the same p(t) to the model both in TVB and MATLAB. 

Thank you,

[Marmaduke Woodman]

Le 27 juil. 2016 à 10:59, Toygun Başaklar <tbas...@gmail.com> a écrit :

Sigma is the level of noise, and this model requires per state variable dispersion values. (so, it’s not the input, technically speaking)

But in terms of JR 1995, it is considered as the input to the system, isn't it? ("p(t)" in the below figure is sigma[3] in the code as far as I understood). I'd like to stimulate the model as described in the paper itself. In order to do that, I should add the desired stimulus (flash pulse in this case) to the noise but since the noise is not a data array (it is a noise object), I cannot add the stimulus. Can we somehow obtain access to this noise object's content (time,data) and manipulate it? 

In TVB, the noise enters into the integration scheme in proportion to the square root of the time step, but if you still want to do it, you can look at the source for the Noise classes and generate arbitrary data.

 
Another question that appeared to my mind, is applying a regional stimulus to this model in TVB equivalent to JR's way of applying stimulus (adding exponential onto white noise) in their paper?

Which section of the paper is that? I am looking but don’t see a description of the noise process they used. 

OTOH,  they’re using a RK4-5 scheme, so if you are using the same in MATLAB, this will already be a source of discrepancy with TVB, since this scheme is not implemented. 

[Toygun Başaklar]

In TVB, the noise enters into the integration scheme in proportion to the square root of the time step, but if you still want to do it, you can look at the source for the Noise classes and generate arbitrary data.

Alright, thank you. I really couldn't find out how to add a transient signal on top of "noise = numpy.sqrt(self.dt) * self.random_stream.normal(size=shape)".

I tried stimulating regions 35 and 36 with the following stimulus

The resulting EEG does not seem to respond to the stimulus at all and it makes me think I am doing something wrong. This may be due to two things:

1- What JR observed in their paper is not EEG recording and we can not distinguish the response of the system inside the EEG.

2- This way of applying input does not produce the same results. (*)
 

Which section of the paper is that? I am looking but don’t see a description of the noise process they used. 

OTOH,  they’re using a RK4-5 scheme, so if you are using the same in MATLAB, this will already be a source of discrepancy with TVB, since this scheme is not implemented. 

It is in section 3.3 Evoked Potentials, but yes we are using the same in MATLAB so I guess they (applying a regional stimulus and JR's way of applying stimulus) are not equivalent. However I think this is irrelevant as noise and stimulus enter the scene in the same equation in the integrator and this also proves the above argument (*) irrelevant.

Below is the code I used if you want to take a look:

from pylab import *
from tvb.simulator.lab import *
import numpy
import tvb.datatypes.projections as projections

jrm = models.JansenRit(v0=6.0)

phi_n_scaling = (jrm.a * jrm.A * (jrm.p_max-jrm.p_min) * 0.5 )**2 / 2.
sigma         = numpy.zeros(6) 
sigma[3]      =

 phi_n_scaling

conn = connectivity.Connectivity(load_default=True)

weighting = numpy.zeros((76, ))
weighting[[35, 36]] = 3.0

eqn_t = equations.PulseTrain()
eqn_t.parameters["onset"] = 5e3 # ms
eqn_t.parameters["tau"]   = 100.0   # ms
eqn_t.parameters["T"]     = 5e3  #
eqn_t.parameters["amp"]   = 200.0

stim = patterns.StimuliRegion(temporal = eqn_t,
                                  connectivity = conn,
                                  weight = weighting)
stim.configure_space()
stim.configure_time(numpy.arange(0., 1e4, 2**-4))

plot_pattern(stim)
ion()
show()

pr = projections.ProjectionSurfaceEEG(load_default=True)
ss = sensors.SensorsEEG.from_file(source_file="eeg_brainstorm_65.txt")
rm = region_mapping.RegionMapping(load_default=True)

rec = (monitors.EEG(projection=pr, sensors=ss, region_mapping=rm, period=1e3/2048.))

sim = simulator.Simulator(
    model=jrm,
    connectivity=conn,
    coupling=coupling.SigmoidalJansenRit(),
    integrator=integrators.HeunStochastic(dt=2 ** -4, noise=noise.Additive(nsig=sigma)),
    monitors= rec,
    stimulus = stim,
    simulation_length=1e4,).configure()

(time, data), = sim.run()

figure()
plot(time, data[:, 0, 31, 0], 'k', alpha=1)
title("EEG")
ioff()
show()

[Marmaduke Woodman]

Le 27 juil. 2016 à 15:55, Toygun Başaklar <tbas...@gmail.com> a écrit :

It is in section 3.3 Evoked Potentials,

OK, sorry I hadn’t read this paper fully before. 

Indeed, their stimulation paradigm is more complex (in particular the randomly distributed impulses) than what we have in our tutorials, but you should be able to accomplish this by implementing your own Stimulus subclass. The stimulus is called by the simulator to produce a stimulation pattern on each time step, 

    stimulus[self.model.cvar, :, :] = self.stimulus(stim_step).reshape((1, -1, 1))

On each time step, evaluate both the random impulse and transient for the regions you want to stimulate and return a suitable vector. 

If you need some help with this, let us know.

Cheers,

Marmaduke 

[Toygun Başaklar]

OK, sorry I hadn’t read this paper fully before. 

Indeed, their stimulation paradigm is more complex (in particular the randomly distributed impulses) than what we have in our tutorials, but you should be able to accomplish this by implementing your own Stimulus subclass. The stimulus is called by the simulator to produce a stimulation pattern on each time step, 

    stimulus[self.model.cvar, :, :] = self.stimulus(stim_step).reshape((1, -1, 1))

On each time step, evaluate both the random impulse and transient for the regions you want to stimulate and return a suitable vector. 

If you need some help with this, let us know.

I would like to clarify/ask a few things: (please correct me if I am wrong)

1- I need to use regional stimulus and stimulate occipital areas as I did before.
2- I will create the transient as the temporal pattern of the stimulus (an arbitrary eqn_t in this case).
* I can add the noise onto eqn_t as well. For this I need to create an arbitrary TemporalApplicableEquation subclass in equations.py.
* I can add the random noise onto this stimulus in the simulator, at the line you provided in the above post (in _loop_update_stimulus).

Which one of the above statements are more suitable?

3- Also, I will not use per state variable dispersion, as I already will include noise in my stimulus. This leads to another point, which integrator scheme should I use?

Please do warn me if there is anything else I should consider.

[Marmaduke Woodman]

Le 28 juil. 2016 à 10:23, Toygun Başaklar <tbas...@gmail.com> a écrit :

I would like to clarify/ask a few things: (please correct me if I am wrong)

1- I need to use regional stimulus and stimulate occipital areas as I did before.
2- I will create the transient as the temporal pattern of the stimulus (an arbitrary eqn_t in this case). 
* I can add the noise onto eqn_t as well. For this I need to create an arbitrary TemporalApplicableEquation subclass in equations.py.
* I can add the random noise onto this stimulus in the simulator, at the line you provided in the above post (in _loop_update_stimulus).

Which one of the above statements are more suitable?

Maybe I need another coffee, but I didn’t see the conflict among those statements; I think they are all correct. 

3- Also, I will not use per state variable dispersion, as I already will include noise in my stimulus. This leads to another point, which integrator scheme should I use?

Heun or Euler deterministic scheme, depending on what you’re using in MATLAB.

[Toygun Başaklar]

28 Temmuz 2016 Perşembe 11:52:57 UTC+3 tarihinde mmwoodman yazdı:

Le 28 juil. 2016 à 10:23, Toygun Başaklar <tbas...@gmail.com> a écrit :

I would like to clarify/ask a few things: (please correct me if I am wrong)

1- I need to use regional stimulus and stimulate occipital areas as I did before.
2- I will create the transient as the temporal pattern of the stimulus (an arbitrary eqn_t in this case). 
* I can add the noise onto eqn_t as well. For this I need to create an arbitrary TemporalApplicableEquation subclass in equations.py.
* I can add the random noise onto this stimulus in the simulator, at the line you provided in the above post (in _loop_update_stimulus).

Which one of the above statements are more suitable?

Maybe I need another coffee, but I didn’t see the conflict among those statements; I think they are all correct.

Let me better explain what I meant in those two statements:

In the 1st statement, say, I have a regional stimulus with temporal equation eqn_t. I can manipulate this eqn_t to introduce noise to my stimulus.

In the 2nd statement, again say, I have a regional stimulus with temporal equation eqn_t. I can manipulate the function _loop_update_stimulus() of the simulator class to introduce noise at each time step.

I wonder which method will be more feasible/meaningful to implement?

Heun or Euler deterministic scheme, depending on what you’re using in MATLAB.

Will use Euler deterministic scheme.

Thank you once more for your time.

[Marmaduke Woodman]

Le 28 juil. 2016 à 12:43, Toygun Başaklar <tbas...@gmail.com> a écrit :

In the 1st statement, say, I have a regional stimulus with temporal equation eqn_t. I can manipulate this eqn_t to introduce noise to my stimulus.

In the 2nd statement, again say, I have a regional stimulus with temporal equation eqn_t. I can manipulate the function _loop_update_stimulus() of the simulator class to introduce noise at each time step.

I wonder which method will be more feasible/meaningful to implement?

Ah I see. Definitely don’t modify _loop_update_stimulus(), instead use approach (1), provide a stimulus class which returns the desired values.

[Toygun Başaklar]

Hello again,

Ah I see. Definitely don’t modify _loop_update_stimulus(), instead use approach (1), provide a stimulus class which returns the desired values.

I've been working on this for some time now. I managed to create a stimulus subclass and apply the same stimulus as in matlab.


Left image is the TVB input, middle uppermost is the MATLAB input. However, TVB does not seem to respond to the stimulus at all. Middle lower plots are MATLAB output, rightmost image is the TVB output.
Below is the code I used.

from pylab import *
from tvb.simulator.lab import *
import numpy
import tvb.datatypes.projections as projections

jrm = models.JansenRit(v0=6.0)

conn = connectivity.Connectivity(load_default=True)

weighting = numpy.zeros((76, ))
weighting[[35,36]] = 3.5

eqn_t = equations.JansenRitExp()
eqn_t.parameters['onset'] = 5000

stim = patterns.StimuliRegion(temporal = eqn_t,
                                  connectivity = conn,
                                  weight = weighting)

pr = projections.ProjectionSurfaceEEG(load_default=True)

ss = sensors.SensorsEEG.from_file(source_file="eeg_brainstorm_65.txt")
rm = region_mapping.RegionMapping(load_default=True)

rec = (monitors.EEG(projection=pr, sensors=ss, region_mapping=rm, period=1))

sim = simulator.Simulator(
    model=jrm,
    connectivity=conn,
    coupling=coupling.SigmoidalJansenRit(),

    integrator=integrators.EulerDeterministic(dt=1),

    monitors= rec,
    stimulus = stim,

    simulation_length=1e4).configure()

(time, data), = sim.run()

plot_pattern(sim.stimulus)
ion()
show()


figure()
plot(time, data[:, 0, 60, 0], 'k', alpha=1)

title("EEG")
ioff()
show()

I am not sure which values to use for dt and monitor period. When I changed from seconds to miliseconds in MATLAB (dt was in seconds in my MATLAB implementation), the model started to behave incorrectly (but it behaved correctly when dt is set to dt/1000 in the integration loop). Is there any way to define dt and time arrays in terms of seconds in TVB?

Also, what we plot in MATLAB is y1-y2 (y_n being the state variables). Can we plot the same in TVB? Do we have access to state variables directly? How is the EEG signal we plot in TVB related to state variables?

Another thing is there are variables in jansen_rit.py named as pmax, pmin and mu. Are these functional? These seem to be used to generate the random input JR talks about.

One last thing, I do not understand the data[:, 0, 60, 0] notation. What do 0, 60 and 0 stand for here? I guessed 60 is for channel selection, but how do I know what channel corresponds to 60?

Is there a flowchart related to how TVB simulator works? I mean I sometimes get confused while tracking the code from one class to another. A generic simulator flowchart (what are initiated inside and with which order, when is integrator used etc.) would be useful.

Thank you.


[Marmaduke Woodman]

Le 1 août 2016 à 10:14, Toygun Başaklar <tbas...@gmail.com> a écrit :

 Is there any way to define dt and time arrays in terms of seconds in TVB?

No, everything is in milliseconds in the simulator. This is a frequent source of issues.

Also, what we plot in MATLAB is y1-y2 (y_n being the state variables). Can we plot the same in TVB? Do we have access to state variables directly? How is the EEG signal we plot in TVB related to state variables?

You can access the state variables directly using the Raw monitor, in which case you can extract y1 and y2 directly. 

The EEG monitor should produce certain variables of interest, of which y1 - y2 should be one, for the Jansen Rit model, but current isn’t the case. By default, EEG data for the JR model should have y0 to y3, and given the linearity of the EEG forward model, you can produce y1 - y2 quantity yourself with the code data[:, 1] - data[:, 2].

Alternatively, when you create the model, you can customize the variables of interest which are used, e.g.

jr = JansenRit(…, variables_of_interest=[‘y1 - y2’])

in which case, the EEG data should contain only one set of signals, corresponding to this quantity.

Another thing is there are variables in jansen_rit.py named as pmax, pmin and mu. Are these functional? These seem to be used to generate the random input JR talks about.

p_min and p_max are not used in the model definition (the dfun method). It is likely that the parameters were listed directly from the paper before realizing that these in particular would be handled by other components of the TVB simulator. So, no, they are not currently operational.

(I didn’t implement JR in TVB, hence my ignorance here.)

One last thing, I do not understand the data[:, 0, 60, 0] notation. What do 0, 60 and 0 stand for here? I guessed 60 is for channel selection, but how do I know what channel corresponds to 60?



TVB’s time series are always 4 dimensional, where the dimensions are time, variable, node or sensor, mode.

Here, the first zero is to select the first variable of interest, the second zero selects the first mode (JR doesn’t have modes so this must be zero). 60 selects the 60th EEG sensor. 

The information about sensors is stored in a sensors attribute of the corresponding monitor, so take  look there to see which sensor it is.

If you’ve chosen a sensor not near the site of stimulation, that may explain the lack of response.

Is there a flowchart related to how TVB simulator works? I mean I sometimes get confused while tracking the code from one class to another. A generic simulator flowchart (what are initiated inside and with which order, when is integrator used etc.) would be useful.

This is a just request; unfortunately, no such diagram exists. Doubly unfortunate, the code involving stimuli is especially convoluted. When in doubt, try to read the code and don’t hesitate to send us emails. 

Cheers,

Marmaduke

[Toygun Başaklar]

Thank you for the detailed answer. I am now plotting y1-y2 and scaled the stimulus to 0.12 - 0.32 range and now getting better (better meaning more similar to MATLAB) results. I took a look at eeg_brainstorm_65.txt and now I'm using 32nd node (Oz) from which I recorded most consistent data (experimentally). However I still have some questions.

What is the difference between using an all-zero weighting matrix in stimulus and not using/giving any stimulus at all? Both of them still produces oscillations (I also would like to learn why/how?), but these oscillations are vastly different.

Can we use different models in different nodes in the same simulator? Why I want to know is because I would like to use two nodes (one in occipital region the other in prefrontal region as JR discussed) and include the connectivity between these two nodes ONLY (is this doable/possible?). At these two nodes JR model will be used but with different C (J in TVB) parameters (135 in occ. 108 in prefrontal). This will be similar to what JR explained in their paper and would also introduce "intercolumn connectivity constants" (K_1 and K_2 in JR-1995).

[Marmaduke Woodman]

Le 1 août 2016 à 16:48, Toygun Başaklar <tbas...@gmail.com> a écrit :

What is the difference between using an all-zero weighting matrix in stimulus and not using/giving any stimulus at all? Both of them still produces oscillations (I also would like to learn why/how?), but these oscillations are vastly different.

I’m not sure I understand entirely but send some code so I see what you mean.

Can we use different models in different nodes in the same simulator? Why I want to know is because I would like to use two nodes (one in occipital region the other in prefrontal region as JR discussed) and include the connectivity between these two nodes ONLY (is this doable/possible?). At these two nodes JR model will be used but with different C (J in TVB) parameters (135 in occ. 108 in prefrontal). This will be similar to what JR explained in their paper and would also introduce "intercolumn connectivity constants" (K_1 and K_2 in JR-1995).

If it is just the parameter values which vary spatially, yes, you just need to make a vector instead of scalar for the parameter. 

If you want a different set of equations on different nodes, this is not current supported though we are brainstorming how to do it in future versions.

Don’t hesitate if you have other questions,

Marmaduke

[Toygun Başaklar]

I’m not sure I understand entirely but send some code so I see what you mean.

I am asking the difference between this (all-zero weighting) (see bolded lines):

weighting = numpy.zeros((76, ))

eqn_t = equations.JansenRitExp()
eqn_t.parameters['onset'] = 5000

eqn_t.parameters['pulse'] = 1

stim = patterns.StimuliRegion(temporal = eqn_t,
                                  connectivity = conn,

                                  weight = weighting).
.
.
.

sim = simulator.Simulator(
    model=jrm,
    connectivity=conn,
    coupling=coupling.SigmoidalJansenRit(),
    integrator=integrators.EulerDeterministic(dt=1),
    monitors= rec,
    stimulus = stim,
    simulation_length=1e4).configure()

and this (no stimulus):

sim = simulator.Simulator(
    model=jrm,
    connectivity=conn,
    coupling=coupling.SigmoidalJansenRit(),
    integrator=integrators.EulerDeterministic(dt=1),
    monitors= rec,

    stimulus = None,
    simulation_length=1e4).configure()

The first one generates:

 The second one generates:

In MATLAB, when I provide no input at all (no noise or pulse density) I do not see any oscillations. Why is this not the case in TVB?

If it is just the parameter values which vary spatially, yes, you just need to make a vector instead of scalar for the parameter. 

This is a part of what I want, yes. I want two different values for J, a and b running in two different JR models on two nodes only. How can I achieve this?
I'm not sure how TVB connectivity really works, but I'm assuming every node is connected with each other with the coupling equation we choose. I would like to have two active nodes (so if there are 76 nodes in TVB, I want only 2 connected nodes and the remaining 74 unconnected/infinite distance) and on these nodes I want the same set of equations but with different parameters. I hope I could make myself clear, I am having a hard time to express myself - I am overwhelmed by the features/functions in TVB.

Thank you.

 

[Marmaduke Woodman]

The first one generates:<Auto Generated Inline Image 1.png>

 The second one generates:<Auto Generated Inline Image 2.png>

In MATLAB, when I provide no input at all (no noise or pulse density) I do not see any oscillations. Why is this not the case in TVB?

ok, I see. In terms of stimulus, no difference. I would guess the difference between the figures is due to random initial conditions.

It is likely there are oscillations simply because there’s a default non-zero coupling strength.

If it is just the parameter values which vary spatially, yes, you just need to make a vector instead of scalar for the parameter. 

This is a part of what I want, yes. I want two different values for J, a and b running in two different JR models on two nodes only. How can I achieve this?
I'm not sure how TVB connectivity really works, but I'm assuming every node is connected with each other with the coupling equation we choose. I would like to have two active nodes (so if there are 76 nodes in TVB, I want only 2 connected nodes and the remaining 74 unconnected/infinite distance) and on these nodes I want the same set of equations but with different parameters. I hope I could make myself clear, I am having a hard time to express myself - I am overwhelmed by the features/functions in TVB.

OK, then maybe you just want to start with a 2 node network? For example, 

conn = connectivity.Connectivity()

conn.weights = numpy.random.randn(2, 2)

conn.delays = numpy.zeros((2, 2))

jr = JansenRit()

jr.J = numpy.r_[J_val_node_1, J_val_node_2]

# etc for jr.a, jr.b

sim = Simulator(model=jr, connectivity=conn, …)

and then configure your simulator as usual. Does this make more sense?

Marmaduke

[Toygun Başaklar]

OK, then maybe you just want to start with a 2 node network? For example, 

conn = connectivity.Connectivity()

conn.weights = numpy.random.randn(2, 2)

conn.delays = numpy.zeros((2, 2))

jr = JansenRit()

jr.J = numpy.r_[J_val_node_1, J_val_node_2]

# etc for jr.a, jr.b

sim = Simulator(model=jr, connectivity=conn, …)

and then configure your simulator as usual. Does this make more sense?

Yes, this is how I guessed it would be done. Thank you.
 
So, K1 and K2 in JR1995 are now the antidiagonal elements of the weight matrix, right? In MATLAB, say we use 4000 and 100 respectively for these values.
To achieve the same, do we set conn.weights to [[0, 4000],[100,0]]?
 
Also, JR uses 30s^-1 as delay between nodes. Do I set conn.delays to [[0, 0.03],[0.03, 0]] assuming the delays are all in ms in TVB?

Thanks.

[Marmaduke Woodman]

Le 2 août 2016 à 09:39, Toygun Başaklar <tbas...@gmail.com> a écrit :

To achieve the same, do we set conn.weights to [[0, 4000],[100,0]]?
 

yes, but make sure it’s an array and not nested lists 

conn.weights = np.array([[0.0, 4000.0, [100.0, 0.0]])

Also, JR uses 30s^-1 as delay between nodes. Do I set conn.delays to [[0, 0.03],[0.03, 0]] assuming the delays are all in ms in TVB?

By 30s^-1, you mean a delay of 1/30 s? If so, this would be

conn.delays = np.array([[0.0, 1e3/30.0], [1e3/30.0, 0.0]])

[Toygun Başaklar]

2 Ağustos 2016 Salı 11:34:32 UTC+3 tarihinde mmwoodman yazdı:

I am using JR's notation with 30s^-1. It's just before the Results section in JR1995 if you want to take a look.

We have started with:

conn = connectivity.Connectivity()
conn.weights = numpy.r_['0,2',[0.0, 4000.0],[100.0,0.0]]
conn.delays = numpy.r_['0,2',[0.0, 0.03],[0.03,0.0]]

and received errors related to centres and orientations. Then we chose one node from prefrontal cortex which is 19th node and one from the visual cortex which is 35th node and modified the code like below;

conn.weights = numpy.r_['0,2',[0.0, 4000.0],[100.0,0.0]]
conn.delays = numpy.r_['0,2',[0.0, 0.03],[0.03,0.0]]
conn.orientations = numpy.r_['0,2',[-0.17818834,  0.66607632,  0.00588274],[-0.09655664,  0.15027455, -0.04954971]]
conn.centres = numpy.r_['0,2',[-70.449382, -12.870265, -11.423195],[ 59.612735,  -2.407101,  -0.844581]]

and received the following error

  File "/home/ee213/Desktop/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb/simulator/simulator.py", line 304, in configure
    self._configure_monitors()
  File "/home/ee213/Desktop/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb/simulator/simulator.py", line 538, in _configure_monitors
    monitor.config_for_sim(self)
  File "/home/ee213/Desktop/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb/simulator/monitors.py", line 586, in config_for_sim
    super(EEG, self).config_for_sim(simulator)
  File "/home/ee213/Desktop/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb/simulator/monitors.py", line 481, in config_for_sim
    numpy_add_at(gain.T, self.rmap, self.gain.T)
IndexError: index 36 is out of bounds for axis 0 with size 2

Could you help me to resolve this issue? Are we on the right track?

Thank you.

[Marmaduke Woodman]

Le 3 août 2016 à 10:15, Toygun Başaklar <tbas...@gmail.com> a écrit :

Could you help me to resolve this issue? Are we on the right track?

The error is due to the precomputed EEG gain matrix expecting 76 regions, while you have but 2; either the precomputed matrix needs to be reduced or you can use the analytic gain matrix.

I think you can load the default EEG sensors, and build the EEG monitor like so:

eeg_sens = sensors.EEG(load_default=True)

eeg_mon = monitors.EEG(sensors=eeg_sens)

and it should calculate an analytic EEG gain matrix (following Sarvas 87).

Can you try that? I don’t have a TVB setup handy to test but if it doesn’t work, I’ll work up a working example.

cheers,

Marmaduke

[Toygun Başaklar]

Without region mapping it didn't work. But with following it worked:

rec = (monitors.EEG(sensors=sensors.SensorsEEG(load_default=True), region_mapping=rm, period=1))

So I omitted the projections object from the monitor. Is there any chance it would cause any problem/error?

[Marmaduke Woodman]

Le 3 août 2016 à 11:55, Toygun Başaklar <tbas...@gmail.com> a écrit :

Without region mapping it didn't work. But with following it worked:

rec = (monitors.EEG(sensors=sensors.SensorsEEG(load_default=True), region_mapping=rm, period=1))

Great! I knew I was forgetting something..

So I omitted the projections object from the monitor. Is there any chance it would cause any problem/error?

No that’s what I meant; in the absence of a provided gain matrix, it will compute a new one for the source/sensor combinations, with an analytic formula.

--
You received this message because you are subscribed to the Google Groups "TVB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tvb-users+...@googlegroups.com.

To post to this group, send email to tvb-...@googlegroups.com.
Visit this group at https://groups.google.com/group/tvb-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/9743ac19-ef09-4903-b1b4-21b95faee87d%40googlegroups.com.

[Toygun Başaklar]

Oh I see. Thank you.

Now, I want to plot the outputs (y1-y2) of both nodes to verify I obtain the same results as JR1995(I will sweep K1 and K2 values). I did the following while instantiating the JR model:

jrm = models.JansenRit(v0=6.0, variables_of_interest=['y1-y2','y1-y2'])

The simulation output now includes two results, but when I plot them

output1 = data[:, 1, 32, 0]
figure('Visual Cortex Node')
plot(time, output1, 'k', alpha=1)
title("EEG")
ion()
show()

output2 = data[:, 0, 32, 0]
figure('Prefrontal Cortex Node')
plot(time, output2, 'k', alpha=1)
title("EEG")
ion()
show()

the plots are identical. Is this because I am recording the same output with the above notation, or could there be something wrong with the simulations?

[Marmaduke Woodman]

32 is the index of the sensor here. Also, you’ve made the same variables of interest, so in this case, you’re plotting the same variable of interest (y1 - y2) for the same sensor. 

You could use a single variable of interest [‘y1 - y2’] and then plot different sensors.

--
You received this message because you are subscribed to the Google Groups "TVB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tvb-users+...@googlegroups.com.
To post to this group, send email to tvb-...@googlegroups.com.
Visit this group at https://groups.google.com/group/tvb-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/34a74999-3ffe-44f9-b510-bb8095863a9b%40googlegroups.com.

[Toygun Başaklar]

32 is the index of the sensor here. Also, you’ve made the same variables of interest, so in this case, you’re plotting the same variable of interest (y1 - y2) for the same sensor. 

You could use a single variable of interest [‘y1 - y2’] and then plot different sensors.

Oh, my bad. I still think as if I'm using MATLAB.

So in this case, I will need to plot from two sensors, one from occipital region and one from prefrontal region, to see the outputs from each column/node. I chose regions 60 (Oz) and 57 (Fpz) and plot the results however the results do not coincide with JR's plots. I have two questions:
1- Do conn.weights correspond to K_1 and K_2 (intercolumn connectivity constants) in JR paper? Should their magnitude be the same as in the paper (K_1 = 4000, K_2 = 100) ? I am asking because I cannot think of any other reason for seeing wrong/irrelevant results (you can see the code used below).
2- Can I give two different stimulus to different nodes? I want to apply noise to prefrontal node and noise+pulse to visual node.

from pylab import *
from tvb.simulator.lab import *
import numpy
import tvb.datatypes.projections as projections

# # Simulate with the Jansen and Rit model.
jrm = models.JansenRit(v0=6.0, variables_of_interest=['y1-y2'])
jrm.J = numpy.r_[108, 135]
jrm.B = numpy.r_[17.6, 22]

conn = connectivity.Connectivity()
conn.weights = numpy.r_['0,2',[0.0, 4000.0],[100.0,0.0]]
conn.delays = numpy.r_['0,2',[0.0, 0.03],[0.03,0.0]]

conn.orientations = numpy.r_['0,2',[ 0.31577083,  0.43790701, -0.34546315],[-0.09655664,  0.15027455, -0.04954971]]
conn.centres = numpy.r_['0,2',[ 59.286775,  -2.547348, -21.924281],[ 59.612735,  -2.407101,  -0.844581]]

weighting = numpy.zeros((2, ))
weighting[[0,1]] = 1

eqn_t = equations.JansenRitExp()
eqn_t.parameters['onset'] = 5000

eqn_t.parameters['pulse'] = 0

stim = patterns.StimuliRegion(temporal = eqn_t,
                                  connectivity = conn,

                                  weight = weighting)


pr = projections.ProjectionSurfaceEEG(load_default=True)
ss = sensors.SensorsEEG.from_file(source_file="eeg_brainstorm_65.txt")
rm = region_mapping.RegionMapping(load_default=True)

rec = (monitors.EEG(sensors=sensors.SensorsEEG(load_default=True), region_mapping=rm, period=1))

# the other aspects of the simulator are standard

sim = simulator.Simulator(
    model=jrm,
    connectivity=conn,
    coupling=coupling.SigmoidalJansenRit(),
    integrator=integrators.EulerDeterministic(dt=1),
    monitors= rec,
    stimulus = stim,
    simulation_length=1e4).configure()

# run it

(time, data), = sim.run()

plot_pattern(sim.stimulus)
ion()
show()

output1 = data[:, :, 60, 0]#Oz

figure('Visual Cortex Node')

plot(time[1000:10000], output1[1000:10000], 'k', alpha=1)

title("EEG")
ion()
show()

output2 = data[:, :, 57, 0]#Fpz

figure('Prefrontal Cortex Node')

plot(time[1000:10000], output2[1000:10000], 'k', alpha=1)

[Marmaduke Woodman]

Le 3 août 2016 à 16:59, Toygun Başaklar <tbas...@gmail.com> a écrit :

So in this case, I will need to plot from two sensors, one from occipital region and one from prefrontal region, to see the outputs from each column/node. I chose regions 60 (Oz) and 57 (Fpz) and plot the results however the results do not coincide with JR's plots. I have two questions: 
1- Do conn.weights correspond to K_1 and K_2 (intercolumn connectivity constants) in JR paper? Should their magnitude be the same as in the paper (K_1 = 4000, K_2 = 100) ? I am asking because I cannot think of any other reason for seeing wrong/irrelevant results (you can see the code used below).

A few guesses:

- The EEG sensors and source nodes are not in the correct positions. One could check the raw time series of the nodes to see if this is the source of discrepancy

- The JR equations as implemented in TVB are different / adapted, so that either parameter scales or configurations also need adapting

From the code I’m not 100% sure where the difference might arise, so I’ll have to take a closer look. If you want to share an excerpt of your MATLAB code and results (by private communication if you prefer), it would accelerate the process.

2- Can I give two different stimulus to different nodes? I want to apply noise to prefrontal node and noise+pulse to visual node.

Not with the current stimuli classes, you’d need a custom one, which knows how to combine these features. I can make an example of this.

cheers,

Marmaduke

[Toygun Başaklar]

From the code I’m not 100% sure where the difference might arise, so I’ll have to take a closer look. If you want to share an excerpt of your MATLAB code and results (by private communication if you prefer), it would accelerate the process.

Sure thing, I will send you my MATLAB script via e-mail. Thank you for your interest in the subject, this is amazing help.

 

Not with the current stimuli classes, you’d need a custom one, which knows how to combine these features. I can make an example of this.

It would be perfect.

[Joe Tharayil]

Hi,

Is there now a publicly available method to do this? I'm trying to solve the same problem.

Thanks,

Joe

[WOODMAN Michael]

hi

Please restate the problem your trying to solve, or quote a more precise section of a previous email.

thanks,

Marmaduke

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Joe Tharayil <thara...@gmail.com>
Sent: Wednesday, August 28, 2019 10:27:40 PM
To: TVB Users
Subject: Re: [TVB] RuntimeWarning in simulate_region_Jansen_Rit

 

--

You received this message because you are subscribed to the Google Groups "TVB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tvb-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/16fb8f35-bfe2-4a4e-ad04-e063c03bec1b%40googlegroups.com.

[Joe Tharayil]

Hi Marmaduke,

I'm trying to replicate the VEP results in the Jansen-Rit paper. If I understand the previous email exchange correctly, the correct way to do this is to define a Stimulus object that applies the background stochastic input to all nodes, and the transient input to the specified nodes, and to run the model with a deterministic integrator. Is this understanding correct? If so, is there code available to define a Stimulus object that applies different stimuli to different nodes, or do I need to write it from scratch?

Best,

Joe

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/9359f128437b46a1a8fbf69edd523468%40univ-amu.fr.

[Gustavo Patow]

Hi,

  I am also trying to reproduce the results in the original JR1995 paper, in particular, Figure 3, which seems to be the definitive test of correctness (and control) of that model. The demo code at

https://nbviewer.jupyter.org/url/docs.thevirtualbrain.org/demos/simulate_region_jansen_rit.ipynb

  is elegant, but does not reproduce that figure... it seemed Toygun Başaklar had it, but no final code posted...

  best

gus.-

PS: I've found a library, called PyRates, that seems to reproduce that figure:

https://github.com/pyrates-neuroscience/PyRates/blob/master/documentation/JansenRit_model_dynamics.ipynb

but I haven't tested it, yet...

To unsubscribe from this group and stop receiving emails from it, send an email to tvb-...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/16fb8f35-bfe2-4a4e-ad04-e063c03bec1b%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "TVB Users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to tvb-...@googlegroups.com.

[Deepa Tilwani]

Hello,

How did you pass ['y1-y2'] ?
To me it shows error.



[WOODMAN Michael]

hi

 Use a list ['y1 - y2'] instead of just a string 'y1-y2'.

cheers,

Marmaduke

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Deepa Tilwani <tilwan...@gmail.com>
Sent: Monday, June 12, 2023 9:50:13 PM
To: TVB Users
Subject: [RESEAUX SOCIAUX] Re: [TVB] RuntimeWarning in simulate_region_Jansen_Rit

 

--
You received this message because you are subscribed to the Google Groups "TVB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tvb-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/95d95bc5-829b-40a6-a3a8-a5ce1a6d053an%40googlegroups.com.