TVB-NEST interface with the JansenRit meanfield model

[Ninontwik]

Hello again,

As mentioned before, I am attempting to include a NEST spiking network model of the thalamus into a TVB whole-brain network which uses the JansenRit meanfield model. I am working on the ParallelNEST31 branch, and using python 3.9.2. I am trying to adapt the minimal Wilson Cowan example notebook such that it will work with the JansenRit TVB model.

The main differences between WilsonCowan and JansenRit are of course the state- and coupling variables. While in Wilson Cowan these variables represent rate, in Jansen Rit we can only work with the post-synaptic potentials (y0, y1 and y2, of which only y1 and y2 are coupling variables). I would like to convert these PSPs to a spike rate that can be inputted into a TVB proxy node (so a NEST spike generating device). As I understand, I can do this by defining a new transformer class, that does this operation from PSP to spikerate. Thus, I am writing this: class Sigmoid(Transformer) in the core.interfaces.base.transformers.models.base.py. And I have also added this new transformer to the DefaultTVBtoSpikeNetTransformers in the builders.py. The operation of the transformer will be very similar to the one in the SigmoidalJansenRit TVB coupling class. Thus, it will need both y1 and y2 as inputs. So in the main file I have defined:

    tvb_spikeNet_model_builder.output_interfaces = \
        [{'voi': np.array(["y1", "y2"]),         # TVB state variables to get data from
          'populations': np.array(["TC"]), # NEST populations to couple to
          'model': 'RATE',                # This can be used to set default tranformer and proxy models
          'coupling_mode': 'TVB',         # or "spikeNet", "NEST", etc
          'proxy_inds': nest_nodes_ids,  # TVB proxy region nodes' indices
          'proxy_model': "RATE",  
          'transformer_model': "SIGMOID",  
          'spiking_proxy_inds': nest_nodes_ids  # Same as "proxy_inds" for this kind of interface
         }
        ]

However, I get this error when trying ot configure the tvb_spikeNet_model_builder:

Exception has occurred: ValueError
2 is not in list
  File "/home/docker/packages/tvb-multiscale/tvb_multiscale/core/interfaces/tvb/interfaces.py", line 60, in <listcomp>
    return np.array([simulator_inds_list.index(ind) for ind in inds])
  File "/home/docker/packages/tvb-multiscale/tvb_multiscale/core/interfaces/tvb/interfaces.py", line 60, in _set_local_indices
    return np.array([simulator_inds_list.index(ind) for ind in inds])
  File "/home/docker/packages/tvb-multiscale/tvb_multiscale/core/interfaces/tvb/interfaces.py", line 67, in set_local_voi_indices
    self.voi_loc = self._set_local_indices(self.voi, monitor_voi)

I think this occurs because self.voi contains the absolute indices of the voi's wrt all the possible vois, so it is [1,2] (y1 and y2), while monitor_voi is set to be the indices of the coupling variables, so [0, 1] because the only coupling variables are y1 and y2.

Am I making a mistake or is something wrong with the code? I can circumvent the error when in core.interfaces.tvb.interfaces.py --> TVBInterface --> set_local_voi_indices, I change self.voi_loc = self._set_local_indices(self.voi, monitor_voi) to 'self.voi_loc = monitor_voi'. But this might cause problems later?

Then the next problem is that I need to use both y1 and y2 in the Sigmoid(Transformer) class to compute the output_buffer. It is unclear to me how to extract y1 and y2 from the input_buffer. If input_buffer would already be (y1-y2), the class would look something like this I think:

class Sigmoid(Transformer):
    """
        Sigmoid Transformer applies a sigmoid function to the input in order to compute the output.
        It comprises of:
            - an input buffer data numpy.array,
            - an output buffer data numpy.array,
            - the parameters of the sigmoid function,
            - a method to apply the sigmoid to the input buffer.

        .. math::
        c_{min} + (c_{max} - c_{min}) / (1.0 + \exp(-a(x-midpoint)/\sigma))
    """
    cmin = NArray(
        label=":math:`c_{min}`",
        default=np.array([0.0,]),
        domain=Range(lo=-1000.0, hi=1000.0, step=10.0),
        doc="Minimum of the sigmoid function",)

    cmax = NArray(
        label=":math:`c_{max}`",
        default=np.array([2.0 * 0.0025,]),
        domain=Range(lo=-1000.0, hi=1000.0, step=10.0),
        doc="Maximum of the sigmoid function",)

    midpoint = NArray(
        label="midpoint",
        default=np.array([6.0,]),
        domain=Range(lo=-1000.0, hi=1000.0, step=10.0),
        doc="Midpoint of the linear portion of the sigmoid",)

    r  = NArray(
        label=r":math:`r`",
        default=np.array([1.0,]),
        domain=Range(lo=0.01, hi=1000.0, step=10.0),
        doc="the steepness of the sigmoidal transformation",)

    a = NArray(
        label=r":math:`a`",
        default=np.array([1.0,]),
        domain=Range(lo=0.01, hi=1000.0, step=10.0),
        doc="Scaling of the coupling term",)

    @property
    def _a(self):
        return self._assert_size("a")

    def configure(self):
        super(Sigmoid, self).configure()
        self._a

    def compute(self):
        """Method that just scales input buffer data to compute the output buffer data."""
        if isinstance(self.input_buffer, np.ndarray):
            self.output_buffer = self.a * self.cmax / (1.0 + np.exp(self.r * (self.midpoint - self.input_buffer)))
        else:
            self.output_buffer = []
            for a, input_buffer in zip(self._a, self.input_buffer):
                self.output_buffer.append(self.a * self.cmax / (1.0 + np.exp(self.r * (self.midpoint - self.input_buffer))))

    def print_str(self):
        return super(Scale, self).print_str() + \
               "\n     - a = %s" % str(self.a)

So the question now is how to obtain y1 and y2 separately to use them properly.

The NESTtoTVB interface will also be a problem since the NEST model will output through a spike recorder, but the TVB node needs to have both it's cvars overwritten which are y1 and y2. And these are difficult to relate to the spiking rate of the pyramidal neurons. One option would maybe be to apply the inverse sigmoid function to this spikerate, and then write this all to y1, and set y2 to zero. Another option would be to just write spikerate to y1 and then change the coupling function of the TVB model specifically for couplings between the thalamus nodes and other nodes.

For the first option, how would I go about writing a transformer that writes specific things to y1 and y2?

I hope the questions are clear and relevant. If I could better ask these questions through another method/platform, please let me know. Also let me know if any additional information could be of help.

Thanks in advance!

[dionperd]

Hello,

Unfortunately, I didn't have time to deal with this problem yet.

I had avoided until now an interface for Jansen-Rit for the complications you are facing now!

Well done by the way to getting that far!

For the TVB -> NEST direction:

I have to see if it is possible to map two state variables from the TVB side, to a single proxy input in the NEST side. Otherwise, it might be better to pass them with separate interfaces and assume one of them to be inhibitory in NEST. Another way is to add a "pseudo" or "non-integrated", intermediate state variable in an augmented version of the Jansen-Rit model, that computes the y1-y2 quantity. Then, you could add that variable to the coupling variables. I have to think more if this would be actually a better solution. Check the ReducedWongWang TVB model for what non-integrated variables are.

As for the vois, it should be functioning. I will have to test it myself to see what is the problem. That method, self._set_local_indices(self.voi, monitor_voi), is supposed to be doing exactly this, i.e., converting an index of state variable that runs for all state variables of a model, to an index that refers only to its coupling variables...

For the NEST -> TVB direction:

Indeed, it is a current limitation that the variables interfacing in the two directions have to be identical, i.e, we cannot avoid updating both of them.

Check how we do this for the Reduced Wong Wang model, where we are even integrating in order to update all state variables of the model.

I find it very interesting to build an interface example for Jansen-Rit. So, I promise to have a look at it once I am back from holiday early next year.

[Ninontwik]

Hello @dionperd ,

I was wondering when you will have a chance to look at the things you've mentioned after your holiday, just so I can adapt my own timeplan a bit. So to see if non-integrated variables are a better solution, if there is something wrong with the vois, and if it is possible to build an interface example for Jansen-Rit.

And please let me know if you need any more information from me.

[dionperd]

Hello. I am dealing with it as for a couple of days already. I might have an example today or tomorrow.

[dionperd]

Hello, have a look at a first example.

On Monday, I will improve it, and provide more explanations.

[Ninontwik]

Thank you so much for this! This is great!

The example is quite clear already but it is a little bit difficult to translate it to my situation.

So just to check if I understand correctly: TVB -> NEST At every TVB node, except for the nodes to be replaced by NEST, a sort of proxy node is created where the Sigmoid function is applied to y1 and y2 to obtain a spikerate that is used as an input to the NEST model?

NEST -> TVB The spikerate of the NEST populations are computed and put through an inverse sigmoid to obtain the PSP's, and with this, y0, y1 and y2 of the TVB nodes are updated. And there is an update from E --> y0, Ein --> y1, and Iin --> y2 ?

I have one main problem with my model. Since I am modeling the Thalamus with NEST, and not a cortical column, it does not make sense to model the neuron populations "E", "Ein", and "Iin". Instead, my NEST model consists of an excitatory population of thalamocortical (TC) neurons, and an inhibitory population of Reticular Nucleus (RN) neurons. The other nodes of the TVB-model project to both these populations. Now, there is a problem when going from NEST to TVB, since these populations do not correspond to y0, y1 and y2 of the Jansen-Rit model. In reality, the TC neurons project to other parts of the brain, the RN neurons receive input from the other parts but project only to the TC neurons. If I'm correct, the coupling of the JansenRit only takes y1 and y2 into account, not y0, so when updating the TVB-thalamus node, the update of y0 is not important, especially not since it can be computed from y1 and y2. Considering these things, I think it would make the most sense to write the spikerate of the "TC" population, through the inverse sigmoid, to y1, and then set y2 to zero.

I haven't figured out how I can best go about doing this. I think the options are:

• make separate interfaces for all connections (something --> y0, "TC" --> y1, and something --> y2) and use the defined "SPIKES" transformer model only for the "TC" --> y1 connections, and write new transformer models for the remaining interfaces, of which the output is just 0.
• make an extra fake neuron population in NEST that is silent (no spiking), and connect this to y0 and y2.

Is the problem clear and do you think one of these options would be reasonable?

Then two other small questions: when defining the parameters of the inverse sigmoid to go from NEST to TVB:

    for interface, N, vv, rr in zip(tvb_spikeNet_model_builder.input_interfaces, 
                                [tvb_spikeNet_model_builder.N_E, 
                                 tvb_spikeNet_model_builder.N_Ein, 
                                 tvb_spikeNet_model_builder.N_Iin], 
                                [-2.5, +20.0, +60.0], 
                                [3.0, 0.75, 0.2]):
        # The "scale_factor" scales the instantaneous rate coming from NEST, before setting it to TVB,
        # in our case converting the rate to a mean reate 
        # and scaling it to be in the TVB model's state variable range [0.0, 1.0]
        interface["transformer_params"] = {"scale_factor": 10*np.array([1e-4])/N,
                                           "Rmax": np.array([12.0]), 
                                           "Rmin": np.array([0.01]),
                                           "v0": np.array([simulator.model.v0[0].item() + vv]),
                                           "r": np.array([simulator.model.r[0].item() * rr])}

What are these values for vv and rr based on? And why are v0 and r computed in this way? I feel like I am missing something about this inverse sigmoid.

And about defining the parameters for the TVB -> NEST interface:

        else:  # RATE
            # Here the rate is a total rate, assuming a number of sending neurons:
            interface["transformer_params"] = {"scale_factor": 100 * np.array([tvb_spikeNet_model_builder.N_Ein])}

Why is the number of neurons multiplied by 100 to obtain the scale factor?

I hope the questions are clear. Please let me know if I can provide more information such as scripts, or if I should provide more information.

[dionperd]

Thank you so much for all the effort! This is great!

I am trying to adapt this code to my own. I get an error when trying: simulator = tvb_spikeNet_model_builder.build() when building the input and output interfaces:

~/packages/tvb-multiscale/tvb_multiscale/core/interfaces/spikeNet/builders.py in _build_tvb_to_spikeNet_interface_proxy_nodes(self, interface)
    238             for i_trg, trg_node in enumerate(interface["spiking_proxy_inds"]):
    239                 weights[i_src, i_trg] = weight_fun(src_node, trg_node, self.tvb_weights)
--> 240                 delays[i_src, i_trg] = delay_fun(src_node, trg_node, self.tvb_delays)
    241                 receptor_type[i_src, i_trg] = receptor_type_fun(src_node, trg_node)
    242                 syn_spec[i_src, i_trg] = syn_spec_fun(src_node, trg_node)

~/packages/tvb-multiscale/tvb_multiscale/core/interfaces/spikeNet/builders.py in _default_tvb_delay_fun(self, source_node, target_node, delays)
    131         if delays is None:
    132             delays = self.tvb_delays
--> 133         return delays[source_node, target_node]
    134 
    135     @abstractmethod

IndexError: index 76 is out of bounds for axis 1 with size 1

Of the TVB model, node 76 and 77 will be replaced by the NEST network. So I think it is trying to get the right delay to or from that node. So source_node = 0 and target_node = 76. However, delays (and also self.tvb_delays) is an array only containing 1 number: [[1.1]]. Even though the TVB simulator: sim.connectivity.delays has a proper matrix with delays connecting all the nodes.

I think the delays are set here in core.interfaces.SpikeNet.builders.py:

    def _get_tvb_delays(self):
        # This is good for ANNarchy because one can set the devices' state at time 0.0
        # For NEST, one has to subtract 1 NEST time step.
        return np.maximum(1,
                          self.tvb_simulator_serialized.get("connectivity.idelays", np.ones((1, 1)))
                          - self.synchronization_n_step + 1) * self.tvb_dt

Which I think will result in one number, namely 2.0. Since (sim_serial is the serialized TVB simulator) sim_serial.get("connectivity.idelays", np.ones((1,1))) gives [[1.1]], synchronizatoin_n_step is 0 and dt = 1.

So I don't know where the [[1.1]] comes from. I think maybe the idelays aren't configured correctly? sim.connectivity.idelays is None, even after configuring the simulator. I cannot find where these idelays should be configured.

Do you have any idea what the problem is here? Is there any other information you need? I could email you my entire code or if you have a hunch I could look into that!

The problem you are describing is due to not having configured your simulator connectivity correctly yet. This is why the serialized simulator doesn't have a tvb_delays entry and sets the default 1D one.

[dionperd]

Thank you so much for this! This is great!

The example is quite clear already but it is a little bit difficult to translate it to my situation.

So just to check if I understand correctly: TVB -> NEST At every TVB node, except for the nodes to be replaced by NEST, a sort of proxy node is created where the Sigmoid function is applied to y1 and y2 to obtain a spikerate that is used as an input to the NEST model?

NEST -> TVB The spikerate of the NEST populations are computed and put through an inverse sigmoid to obtain the PSP's, and with this, y0, y1 and y2 of the TVB nodes are updated. And there is an update from E --> y0, Ein --> y1, and Iin --> y2 ?

I have one main problem with my model. Since I am modeling the Thalamus with NEST, and not a cortical column, it does not make sense to model the neuron populations "E", "Ein", and "Iin". Instead, my NEST model consists of an excitatory population of thalamocortical (TC) neurons, and an inhibitory population of Reticular Nucleus (RN) neurons. The other nodes of the TVB-model project to both these populations. Now, there is a problem when going from NEST to TVB, since these populations do not correspond to y0, y1 and y2 of the Jansen-Rit model. In reality, the TC neurons project to other parts of the brain, the RN neurons receive input from the other parts but project only to the TC neurons. If I'm correct, the coupling of the JansenRit only takes y1 and y2 into account, not y0, so when updating the TVB-thalamus node, the update of y0 is not important, especially not since it can be computed from y1 and y2. Considering these things, I think it would make the most sense to write the spikerate of the "TC" population, through the inverse sigmoid, to y1, and then set y2 to zero.

I haven't figured out how I can best go about doing this. I think the options are:

• make separate interfaces for all connections (something --> y0, "TC" --> y1, and something --> y2) and use the defined "SPIKES" transformer model only for the "TC" --> y1 connections, and write new transformer models for the remaining interfaces, of which the output is just 0.
• make an extra fake neuron population in NEST that is silent (no spiking), and connect this to y0 and y2.

Is the problem clear and do you think one of these options would be reasonable?

Then two other small questions: when defining the parameters of the inverse sigmoid to go from NEST to TVB:

    for interface, N, vv, rr in zip(tvb_spikeNet_model_builder.input_interfaces, 
                                [tvb_spikeNet_model_builder.N_E, 
                                 tvb_spikeNet_model_builder.N_Ein, 
                                 tvb_spikeNet_model_builder.N_Iin], 
                                [-2.5, +20.0, +60.0], 
                                [3.0, 0.75, 0.2]):
        # The "scale_factor" scales the instantaneous rate coming from NEST, before setting it to TVB,
        # in our case converting the rate to a mean reate 
        # and scaling it to be in the TVB model's state variable range [0.0, 1.0]
        interface["transformer_params"] = {"scale_factor": 10*np.array([1e-4])/N,
                                           "Rmax": np.array([12.0]), 
                                           "Rmin": np.array([0.01]),
                                           "v0": np.array([simulator.model.v0[0].item() + vv]),
                                           "r": np.array([simulator.model.r[0].item() * rr])}

What are these values for vv and rr based on? And why are v0 and r computed in this way? I feel like I am missing something about this inverse sigmoid.

And about defining the parameters for the TVB -> NEST interface:

        else:  # RATE
            # Here the rate is a total rate, assuming a number of sending neurons:
            interface["transformer_params"] = {"scale_factor": 100 * np.array([tvb_spikeNet_model_builder.N_Ein])}

Why is the number of neurons multiplied by 100 to obtain the scale factor?

I hope the questions are clear. Please let me know if I can provide more information such as scripts, or if I should provide more information.

I am still in the process of adding more options and documenting the example.

I had to do major changes for this to work.

It would be better if we skyped at some point to explain better the details. I think it would save us some time!

In the mean time, a few elements here:

TVB -> NEST:

Option 1 is when coupling is computed in TVB: tvb_spikeNet_model_builder.default_coupling_mode = TVB_NEST_COUPLING_MODE = 'TVB" # "NEST" # "TVB"

In this case the proper (JansenRit sigmoidal) coupling directed to your NEST nodes is already computed in TVB, with the proper weights and delays determined by the TVB connectome, and the conversion to a rate/current quantity is already made.

So, all you need is a LinearRate transformer to bring it to the correct range of values.

There is 1 proxy TVB node for every NEST spiking node constructed as an inhomogeneus_poisson_generator device that stimulates the target NEST nodes accordingly.

You could also use one of the transformers that produce the spike trains via Elephant.

The trick here is that the interface picks the TVB values from the history buffer to compute the coupling and send it to NEST. JansenRit sigmoidal coupling is particular in that although it takes values y1 and y2, it returns only one value, and therefore we are interested in the position 0 in python indexing of the output coupling. This is why we need to "trick" the interface by asking for 'y0' as the variable of interest of the interface: voi = np.array(["y0"])

Mind also that I had to use all three state variables "y0", "y1", and "y2" as cosimulation vois, because in my example the NEST network is identical to the JansenRit population model, and I want to update all 3 of them in the NEST -> TVB direction. Therefore, the first position is indeed the "y0".

Option 2 is when coupling is computed in NEST. In this case, we need to take y1 and y2 from every node of TVB, compute the JansenRite sigmoidal presynaptic coupling expression as a transformer and then, make 1 TVB proxy node for every TVB node inside NEST, and couple it to the target NEST nodes with weights and delays taken from TVB connectome. This is a slower option since we need to update N-2 proxies, instead of just 2. However, this explains this part of the code: tvb_spikeNet_model_builder._tvb_to_spikeNet_transformer_models = DefaultTVBtoSpikeNetTransformersJansenRitCoupling voi = np.array(["y0", "y1", "y2"]) proxy_inds = np.arange(simulator.connectivity.number_of_regions).astype('i') proxy_inds = np.delete(proxy_inds, nest_nodes_inds)

I have now added such a JansenRit sigmoidal coupling transformer and you can find it in core.interfaces.transformers.models.jansen_rit.py

In both cases the coupling variables of the TVB JansenRit model have to adjusted to return all 3 state variables.

This is now done in a dedicated for tvb-multiscale version of the model: core.tvb.cosimulator.models.jansen_rit.py

In the NEST -> TVB direction there are also now 2 options (pull the latest changes!):

Option 1 is to record spikes, convert them into mean rates and then use an inverse sigmoidal transformer which I have also added in the same core.interfaces.transformers.models.jansen_rit.py

Option 2 is to record voltages directly and average them across the population and scale them accordingly to bring them to the correct range of values of each of y0, y1, y2.

Check in the latest version of the notebook these two options.

The second one returns a more meaningful output but it is much slower because unfortunately we have to return from memory always all the membrane potential data from the multimeters. I will actually try now to record them to files which makes it probably faster to only read the new entries for every synchronization_time, and probably become faster. More on this later! Another difficulty in this case is that to get voltage of time stamp t in NEST you have to simulate also the step t+min_delay.

This is why I start a NEST resolution step before the cosimulation. It is something to discuss with the NEST developers...

Now in your case, it doesn't make sense to substitue the TVB JansenRit data of the two thalamic regions with quantities that have a completely different meaning from your spikiing nodes.

Ideally, you could probably change the mean field model only for those two nodes, even if it is something as simple as:

dx_i/dt = -x_i with x_i being substituted by NEST input for every population i of your spiking thalamic nodes.

This can be done in TVB, i.e., changing the equations for particular nodes with just an IF statement or using np.where, and having a mask across regions.

In case the input is firing rate, you might want to smooth it a bit.

This requires transformers with memory, pretty much like the current Integrator transformer model.

We can further discuss this if necessary.

[dionperd]

A notebook with different kinds of interfaces for JansenRit TVB model can now be found in the tvb_nest examples.