Stimulus on Jansen&Rit

[Rachel S]

Hi,

I am confused by the way stimulus is represented in TVB. For JansenRit model, it is adding the same value to the second and third state variable using stimulus[self.model.cvar, :, :] = self.stimulus(stim_step).reshape((1, -1, 1)). Wouldn't that cause zero influence to the output of the JansenRit, which is y1-y2? Am I missing something? Appreciate your help, thank you.

[Rachel S]

Also, I am confused with the relationship between the stimulus and the noise in the Stochastic Integrator. For many papers, they mentioned p(t) is the input for the model, which can represent the activity from other brain region or external stimuli. I've been playing around with the HeunStochastic Integrator and the stimulus for a while now, and I have a few observations that I cannot understand.

1. In the JansenRit tutorial, integrator noise was 

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

added to y4 (as in the JR equations p(t) is added to y4). 

• What is the idea behind scaling noise with sqrt(dt)? Why not using noise * dt in the integrator.scheme? 
• How is that related to the noise input p(t)?

2. If I use stimulus as p(t) and use deterministic integrator. With a much larger variance, I can get similar behavior as point 1. Why is that? What is the relation between adding noise using Stochastic Integrator and adding noise explicitly by defining the stimuli as Gaussian noise?

Thank you,

Rachel

[WOODMAN Michael]

Hi

 

From the perspective of a neural mass model (especially early papers such as Jansen-Rit), network afferents, noise and stimulus casn often be written as a single “input” term, yet in a modeling project it is interesting to distinguish different mechanisms

 

1. “spontaneous” noise driven fluctuations
2. Coordination through network effects
3. Stimulus driven changes in network activity

 

The TVB software reflects this set of theoretical distinctions, even if you can in practice ignore them for some purposes.

 

Scaling noise by sqrt(dt) is a basic element of SDE theory, please see for example the attached intro.

 

Cheers,

 

Marmaduke Woodman IR AMU INS U1106 +33 7 67 77 84 72

--
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/8588de9d-f5d4-4e79-93cc-0712e2b06e11%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Rachel S]

Hi, 

Thank you for your prompt response. Do you have any comments on my first question about the stimulus for JansenRit?

" For JansenRit model, it is adding the same value to the second and third state variable using stimulus[self.model.cvar, :, :] = self.stimulus(stim_step).reshape((1, -1, 1)). Wouldn't that cause zero influence to the output of the JansenRit, which is y1-y2? Am I missing something? Appreciate your help, thank you. "

Best,

Rachel

WOODMAN Michael <marmaduk...@univ-amu.fr> 于2019年5月7日周二 上午3:30写道：

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/66LFeV0Xtd0/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/e1e025299ff6434fb57b6e95f354c7b5%40univ-amu.fr.

[WOODMAN Michael]

hi

 

I deliberately did not respond in the hopes that someone more knowledgable in Jansen-Rit model would reply...

 

I can speculate that while a stimulation on those pair of variables may be harder to observe, the stimulation would certainly create a perturbation in state space and therefore exert some influence. If you aren't currently observing any effect of the stimulus, it may be worth plotting phase planes of the JR state variables to see what's happening (and refer to the original equations to check magnitidues etc).

 

cheers,

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/CAAx9jGfF9fPLQNb15aT9vaasyXc3xNuw4WSSNX6C4QJOvk9Rjg%40mail.gmail.com.

[Julie Courtiol]

Hi,

I already verified how the stimulus was constructed in TVB. If I remember well, the stimulus is applied only on the first cvar state-variable, here it corresponds to y_1. This has to be re-check again.

Best,

Julie

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

[Gustavo Patow]

Hi,

  Well, actually, the _loop_update_stimulus method from the simulator class reads, as Rachel S pointed:

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

  and the JansenRit class defines that cvar is

 cvar = numpy.array([1, 2], dtype=numpy.int32)

  Which means that the noise is added, through the stimulus construct, to y1 and y2. As Rachel S said, this would give the same stimulus to two vars that are latter subtracted to get the final value we can measure (y1 - y2, if I remember correctly). But I think things are a bit worse, because, in the JansenRit model, any external input should enter trough y4, in the form 

   A a p(t) + other terms

  I guess I am missing something important, as this is NOT what I understand is done through the integrator classes (e.g., the Euler-Maruyama does X = X + dt (dfun + stimulus) + noise). Please, any input of what I am missing here would be greatly appreciated!

  all the best

gus.-

[Joe Tharayil]

Hi Gustavo,

Have you figured this out yet? I've been trying to replicate the original Jansen-Rit paper with TVB, and I haven't been able to do so for the reasons you describe.

Best,

Joe

[WOODMAN Michael]

hi

There are several results in that paper, which are you trying to reproduce?

cheers,

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Joe Tharayil <thara...@gmail.com>
Sent: Tuesday, January 14, 2020 1:56:01 PM
To: TVB Users

Subject: Re: [TVB] Re: Stimulus on 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/ac4e8ed8-853f-41d0-a7a0-585f0bd3da38%40googlegroups.com.

[Joe Tharayil]

Hi,

I'd like to be able to reproduce Figures 3-6, or at least Figures 3 and 5.

Best,

Joe

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

[Gustavo Patow]

Hi Joe,

   I am sorry, at this moment I am not actually using TVB, but my own code, and I was using TVB as a reference. While studying the code I discovered that discrepancy and I decided to report it, mostly to help the community. Meanwhile, I fixed my own code and now I can reproduce both Jansen-Rit 1995 Fig 3 and David-Friston 2003, Fig 4. If you want my (horribly awful) code, I can send you the link to GitHub...

  Note: A colleague of mine is heavily insisting me to move to TVB, so probably I'll end up being a full user soon... ;-)

  all the best

gus.-

[WOODMAN Michael]

Hi

Reading through the paper, the TVB implementation of the single column model is correct, but implementing the two-column case is not easy, because they have added a PSP blocks for the connections between the columns.  In other words, while a single column required 6 ODEs, the two columsn require 2*6 + 2*2 = 16 ODEs.  The straightforward adaptation to TVB would add 2 ODEs per node for the afferent PSP, no longer requiring a difference in the coupling function, and a stochastic integartor should be usable for the noise.

This appears to be done in the contributed Jansen-Rit implementation for the David 2005 paper,

https://github.com/the-virtual-brain/tvb-root/blob/master/scientific_library/contrib/simulator/models/jansen_rit_david.py

but it should be checked again for correctness as it uses another naming scheme so code isn't comparable.

cheers,

Marmaduke

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Gustavo Patow <dag...@gmail.com>
Sent: Wednesday, January 15, 2020 10:13:09 AM
To: TVB Users
Subject: [TVB] Re: Stimulus on 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/c575fd5d-824f-49b6-980e-1dbbb914f9ba%40googlegroups.com.

[WOODMAN Michael]

hi

I've put together the following script,

https://gist.github.com/maedoc/d6780270b47f956c541eef5ed85e8194

which runs a two column simulation with equations as in the paper, by adding an PSP to each node and using its output as the coupling variable.  This is scaled by the connectivity (K1 and K2 here) and then used as input to the JR equation 6.  Furthermore the noise is factored out of the y4 equation by replacing with a mean (mu=0.22) and adding small white noise.

The single column results Fig 3 can be reproduced by changed J (what JR calls C) in the script on line 15, whereas k1 and k2 are below that.  I suspect there's a scaling issue since using k1 & k2 values from the paper don't allow for reproducing all the results shown, but it could be a useful starting point.

cheers,

Marmaduke

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Joe Tharayil <thara...@gmail.com>

Sent: Wednesday, January 15, 2020 10:04:39 AM

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/4297ae98-6820-4c0e-88ad-e01612038ae8%40googlegroups.com.

[Joe Tharayil]

Hi Marmaduke,

Thanks so much for your help. I need to spend a bit of time looking over your script, but one immediate question I have is why it's valid to factor out the input noise in y4_dot, and replace it with white noise in y4. It does seem to work, but shouldn't the input to y4 be Brownian noise instead, since you're effectively integrating over Gaussian noise?

Thanks again!

Best,

Joe

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

[WOODMAN Michael]

Hi

It was just a quick guess that the noise could be pulled out of the RHS of the ODE and applied as a more standard SDE.  After some experimenting, I don't think it works quite well, and I also just did what JR describe by setting all nsig to 0.0 and putting

dy[:6] = super().dfun(state_variables[:6], coupling + np.random.uniform(0.120, 0.320, size=(1, 2, 1)), local_coupling)
       

in the dfun instead, and that seems to be better behaved.  That also lets you use the RK4 in TVB though I didn't notice differences with a dt=0.1 ms (since the dynamics are on the timescale of 100ms, this make sense).

In any case, I'd be happy to see any corrections about the noise term.

cheers,

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Joe Tharayil <thara...@gmail.com>

Sent: Wednesday, January 15, 2020 2:24:32 PM
To: tvb-...@googlegroups.com

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/CADtj9xonhvk4nYHYZhLH-LC-7XX7rq-ssAGU_pMfAHW816VZ%3DQ%40mail.gmail.com.

[WOODMAN Michael]

Hi Rachel

(1) Indeed the papers use p(t) as the input, and this is accomodated in TVB.  What TVB doesn't work so well for is noise terms which are embedded in the RHS equations of the model.  For JR it's

y4' = A a (p(t) + ....)

and the A*a*p(t) can perhaps be factored out and used to scale the noise (with a Gaussian approximation of the uniformly distributed noise) for a proper SDE integration scheme such as HeunStochastic.   This is what is happening in the tutorial, though I think it should be assigned to sigma[4] not sigma[3], since the JR paper counts from y0, not y1.

> What is the idea behind scaling noise with sqrt(dt)?

The noise scales with sqrt(dt) because the equivalent Brownian process' variances increases with the sqrt of time.  Since we are increasing time by dt at each step, the variance is scaled by sqrt(dt). 

> How is that related to the noise input p(t)?

we prefer to formulate the JR equations as SDE since that is what TVB expects, which is why we are trying to find appropriate scaling factors to translate p(t) noise to a normal SDE

> 2. If I use stimulus as p(t) and use deterministic integrator. With a much larger variance, I can get similar behavior as point 1. 

Because the stimulus is expressed along side the RHS of the ODEs for the neural mass, the integrator scales it by dt as well.  This means that if you just trying to inject noise, you need to increase the noise amplitude by dt to achieve the same effect as using stochastic integrator.

Cheers,

Marmaduke 

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Rachel S <sunsunru...@gmail.com>
Sent: Monday, May 6, 2019 11:12:35 PM
To: TVB Users
Subject: [TVB] Re: Stimulus on 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.

[Joe Tharayil]

Hi Marmaduke,

After thinking about this a bit more carefully, I'm really not sure why the output of the additional PSP blocks isn't fed into a sigmoid function before coupling. Doing that would be should be functionally equivalent to using y0 as the coupling variable with TVB's pre-sigmoid coupling function, given the correct delay matrix.

Given all of this, I'm curious about the motivation behind the SigmoidalJansenRit coupling function. Could you point me to the paper that it comes from?

On Wed, Jan 15, 2020 at 1:52 PM WOODMAN Michael <marmaduk...@univ-amu.fr> wrote:

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

[WOODMAN Michael]

hi

The sigmoidal is part of the PSP block, so it makes sense

I didn't implement that component and not aware of any particular paper to be honest. Still, it would be useful if you wanted an instantaneous coupling instead of through 2 ODEs.  I think the JR paper is a challenge to implement with TVB additionally because the connectivity weights K1 and K2 are applied after the PSP blocks, so generalizing to larger networks leaves several choices to the modeler.

cheers,

Marmaduke

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Joe Tharayil <thara...@gmail.com>

Sent: Wednesday, January 15, 2020 4:22:19 PM
To: tvb-...@googlegroups.com

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/CADtj9xrk_PZ_2MpF78vScHTKL5gAmWNdAT%2B1oisnzYi3rUasaw%40mail.gmail.com.

[Rachel S]

Hi Marmaduke,

Thank you so much for your response. I implemented my own version of the Jansen Rit model with heun_deterministic integrator (taking the "dt" term into account as you mentioned earlier) and I am also able to reproduce the single node Jansen rit paper with the same noise parameters described in Wendling, 2000 paper. 

Regarding the noise std, I wanted to add my own observations and wish it would be helpful to the community:

From the simulate_region_jansen_rit.ipynb, the noise added to the system is:

phi_n_scaling = (jrm.a * jrm.A * (jrm.p_max-jrm.p_min) * 0.5 )**2 / 2. 

(BTW, the current demo still adds noise to y[3] instead of y[4])

The reason to calculate the phi_n_scaling in that way is that:

I almost got the same equation other that sqrt(dt) used in the stochastic integrator to stimulate the Brownian process. I assume the math is not as simple as making two terms to be equal and I could make a mistake there.

Anyway, for those who are interested, hope this could start a discussion.

Best regards,

Rachel

[WOODMAN Michael]

hi Rachel

Thanks this is a helpful way of explaining the difference between what the paper does and what TVB does by default.   I think the original set up can be replicated in TVB but requires a uniform noise source & factoring the noise out of the ODEs, which I didn't have time to do before. In any case this should all be enough for someone else to pick it up if desired.

cheers,

Marmaduke Woodman, TVB Engineer, INS AMU; +33 7 67 77 84 72

---

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Rachel S <sunsunru...@gmail.com>

Sent: Wednesday, March 4, 2020 5:27:32 PM
To: TVB Users
Subject: Re: [TVB] Re: Stimulus on 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/54760f58-372e-4ddb-8b9f-90fe3fc0ac43%40googlegroups.com.