Model implementation and execution with Scikit-NeuroMSI

This tutorial covers the basic pipeline for implementing and running multisensory integration models using Scikit-NeuroMSI. We show how to implement the different models classes currently supported by the package, organized in the mle, bayesian and neural modules. It covers the execution objects and methods that are shared by all the models of the package: run and NDResult.

Note: In this tutorial we assume that you already have a basic knowledge of numpy for scientific computing.

Implementation of Maximum Likelihood Estimation (MLE) models

To easily implement the model developed by Alais and Burr (2004) you can import the corresponding module and instatiate the AlaisBurr2004 class:

from skneuromsi.mle import AlaisBurr2004

model = AlaisBurr2004()
model

<skneuromsi.mle._alais_burr2004.AlaisBurr2004 at 0x708e7e48ce50>

Exploration of the model object

By calling the model using the vars function you can explore its main built-in parameters and methods:

vars(model)

{'_mode0': 'auditory',
 '_mode1': 'visual',
 '_position_range': (-20, 20),
 '_position_res': 0.01,
 '_time_range': (1, 1),
 '_time_res': 1.0,
 '_random': Generator(PCG64) at 0x708DC82435A0,
 'run': <function skneuromsi.mle._alais_burr2004.AlaisBurr2004.run(*, auditory_position=-5, visual_position=5, auditory_sigma=3.0, visual_sigma=3.0, noise=None)>}

• mode0: Name of the first sensory modality defined in the model.

• mode1: Name of the second sensory modality defined in the model.

• position_range: The range of possible positions where the stimulus could be delivered (in degrees).

• position_res: The resolution of the range of possible positions where the stimulus could be delivered (in degrees).

• time_range: The range of possible times when the stimulus could be delivered. Here is set to 1 because the model has no temporal dimension.

• time_res: The resolution of the range of possible times when the stimulus could be delivered.

• random: Embedded random number generator (useful to include noise in the model responses).

• run: Executes the model and saves the result.

We can re-implement Alais and Burr (2004) model with different sensory modalities and position ranges by modifying the initial parameters of the model:

model = AlaisBurr2004(
    mode0="visual", mode1="proprioceptive", position_range=(-50, 50)
)
vars(model)

{'_mode0': 'visual',
 '_mode1': 'proprioceptive',
 '_position_range': (-50, 50),
 '_position_res': 0.01,
 '_time_range': (1, 1),
 '_time_res': 1.0,
 '_random': Generator(PCG64) at 0x708DC7FA1540,
 'run': <function skneuromsi.mle._alais_burr2004.AlaisBurr2004.run(*, visual_position=-5, proprioceptive_position=5, visual_sigma=3.0, proprioceptive_sigma=3.0, noise=None)>}

The model object has a built-in run method:

model.run

<function skneuromsi.mle._alais_burr2004.AlaisBurr2004.run(*, visual_position=-5, proprioceptive_position=5, visual_sigma=3.0, proprioceptive_sigma=3.0, noise=None)>

By calling this method we can observe its arguments:

• visual_position: The position where the visual stimulus is delivered (in degrees).

• proprioceptive_position: The position where the propioceptive stimulus is delivered (in degrees).

• visual_sigma: Standard deviation of the auditory estimate.

• proprioceptive_sigma: Standard deviation of the propioceptive estimate.

Now let’s run the model for equidistant auditory and visual locations:

res = model.run(visual_position=-5, proprioceptive_position=5)
res

<NDResult 'AlaisBurr2004', modes=['visual' 'proprioceptive' 'multi'], times=1, positions=10000, positions_coordinates=1, causes=False>

The model outputs one NDResult object containing the results of both unisensory estimators and the multisensory estimator.

Refer to the API documentation for further information about the parameters to manipulate in this model.

Exploration of the results object

By calling the NDResult object using the vars function you can explore its main built-in parameters and methods:

vars(res)

{'_mname': 'AlaisBurr2004',
 '_mtype': 'MLE',
 '_output_mode': 'multi',
 '_nmap': {'auditory': 'visual',
  'visual': 'proprioceptive',
  'auditory_weight': 'visual_weight',
  'visual_weight': 'proprioceptive_weight'},
 '_time_range': array([1, 1]),
 '_position_range': array([-50,  50]),
 '_time_res': 1.0,
 '_position_res': 0.01,
 '_run_parameters': {'visual_position': -5,
  'proprioceptive_position': 5,
  'visual_sigma': 3.0,
  'proprioceptive_sigma': 3.0,
  'noise': None},
 '_extra': {'visual_weight': 0.5, 'proprioceptive_weight': 0.5},
 '_causes': None,
 '_nddata': <xarray.DataArray 'values' (modes: 3, times: 1, positions: 10000,
                             positions_coordinates: 1)> Size: 240kB
 array([[[[1.84356985e-050],
          [1.93808093e-050],
          [2.03741451e-050],
          ...,
          [1.65115301e-074],
          [1.55331368e-074],
          [1.46125559e-074]]],


        [[[1.37463811e-074],
          [1.46125559e-074],
          [1.55331368e-074],
          ...,
          [2.14181549e-050],
          [2.03741451e-050],
          [1.93808093e-050]]],


        [[[4.33458701e-122],
          [4.84392981e-122],
          [5.41300350e-122],
          ...,
          [6.04879857e-122],
          [5.41300350e-122],
          [4.84392981e-122]]]])
 Coordinates:
   * modes                  (modes) <U14 168B 'visual' 'proprioceptive' 'multi'
   * times                  (times) int64 8B 0
   * positions              (positions) int64 80kB 0 1 2 3 ... 9997 9998 9999
   * positions_coordinates  (positions_coordinates) <U2 8B 'x0'}

We observe metadata related to the model that was executed to generate the results object, as well as the model execution (run) parameters. Notably the NDResult object holds useful methods for exploring the results of the model execution: get_modes, stats and plot.

Now let’s explore the get_modes method:

res.get_modes()

		modes	visual	proprioceptive	multi
times	positions	positions_coordinates
0	0	x0	1.843570e-50	1.374638e-74	4.334587e-122
	1	x0	1.938081e-50	1.461256e-74	4.843930e-122
	2	x0	2.037415e-50	1.553314e-74	5.413003e-122
	3	x0	2.141815e-50	1.651153e-74	6.048799e-122
	4	x0	2.251541e-50	1.755135e-74	6.759122e-122
	...	...	...	...	...
	9995	x0	1.865646e-74	2.366862e-50	7.552692e-122
	9996	x0	1.755135e-74	2.251541e-50	6.759122e-122
	9997	x0	1.651153e-74	2.141815e-50	6.048799e-122
	9998	x0	1.553314e-74	2.037415e-50	5.413003e-122
	9999	x0	1.461256e-74	1.938081e-50	4.843930e-122

10000 rows × 3 columns

This method returns a pandas DataFrame with the output values for each modality at each time point and spatial coordinates. In the package we called such output values simply as values, given that the outputs are different across models. In this MLE model, the output values are probability density values.

Note: This model does not include the temporal dimension, hence here is represented as a fixed time point. This model does include spatial dimensions (generically named positions in the package). The package allows to have more than one spatial dimension (in the package called position_coordinates). In this model we have a single spatial dimension, defined as a position_coordinate named x0 by default.

Now let’s explore the stats method:

res.stats()

	describe
count	3.000000e+04
mean	1.000000e-02
std	3.114718e-02
min	4.334587e-122
25%	3.091076e-43
50%	4.385153e-20
75%	3.429870e-06
max	1.880632e-01

By default it provides descriptive statistics of the model output. We can have more specific statistics by specifying the kind and modes argument, for example:

res.stats(kind="max", modes="visual")

0.1329807601338109

A useful method is the dimmax, which provides the maximum value observed in a specific modality and its temporal and spatial coordinates:

res.stats.dimmax(modes="visual")

times                                     0
positions                              4500
positions_coordinates                    x0
values                   0.1329807601338109
Name: max, dtype: object

Refer to the API documentation for further information about the NDResult object.

So far we have covered the fundamentals of model classes and results objects in the Scikit-NeuroMSI package. These are applicable to all model classes. Now we are going to introduce the Bayesian and neural modules, and present some distinctive characteristics of the objects generated by this models.

Implementation of Bayesian models

You can implement the Causal Inference model developed by Kording et al. (2007) by importing the corresponding module and instantiating the Kording2007 class:

from skneuromsi.bayesian import Kording2007

model = Kording2007(position_range=(-20, 20), position_res=1, n=100000)
model

<skneuromsi.bayesian._kording2007.Kording2007 at 0x708db6520bb0>

vars(model)

{'_n': 100000,
 '_mode0': 'auditory',
 '_mode1': 'visual',
 '_position_range': (-20, 20),
 '_position_res': 1.0,
 '_time_range': (0, 1),
 '_time_res': 1.0,
 '_random': Generator(PCG64) at 0x708DC5052C00,
 'run': <function skneuromsi.bayesian._kording2007.Kording2007.run(*, auditory_position=-15, visual_position=15, auditory_sigma=2.0, visual_sigma=10.0, p_common=0.5, prior_sigma=20.0, prior_mu=0, strategy='averaging', noise=True, causes_kind='count', dimension='space')>}

Refer to the API documentation for more details about the available parameters for this model.

Let’s run the model for two conflicting stimulus locations:

res = model.run(auditory_position=-15, visual_position=15, causes_kind="count")
res

<NDResult 'Kording2007', modes=['auditory' 'visual' 'multi'], times=1, positions=40, positions_coordinates=1, causes=2>

In this model we have an output named causes, which provides information about the estimated number of causes (sources) of the stimuli presented to the model. This is a fundamental attribute of causal inference models of multisensory integration.

Let’s now explore the model output:

res.get_modes()

		modes	auditory	multi	visual
times	positions	positions_coordinates
0	0	x0	0.006784	0.001400	0.000000
	1	x0	0.022207	0.005452	0.000012
	2	x0	0.056269	0.018036	0.000061
	3	x0	0.107057	0.051357	0.000230
	4	x0	0.165070	0.102585	0.000460
	5	x0	0.195675	0.163996	0.001102
	6	x0	0.185072	0.200408	0.002312
	7	x0	0.134455	0.189935	0.004043
	8	x0	0.076893	0.139467	0.005690
	9	x0	0.035084	0.078207	0.008062
	10	x0	0.011434	0.034802	0.010011
	11	x0	0.003287	0.011054	0.012614
	12	x0	0.000641	0.002821	0.013558
	13	x0	0.000070	0.000440	0.015180
	14	x0	0.000000	0.000040	0.016972
	15	x0	0.000000	0.000000	0.018618
	16	x0	0.000000	0.000000	0.019018
	17	x0	0.000000	0.000000	0.020591
	18	x0	0.000000	0.000000	0.022480
	19	x0	0.000000	0.000000	0.024017
	20	x0	0.000000	0.000000	0.025083
	21	x0	0.000000	0.000000	0.026632
	22	x0	0.000000	0.000000	0.028375
	23	x0	0.000000	0.000000	0.030615
	24	x0	0.000000	0.000000	0.032612
	25	x0	0.000000	0.000000	0.034634
	26	x0	0.000000	0.000000	0.035905
	27	x0	0.000000	0.000000	0.039827
	28	x0	0.000000	0.000000	0.042272
	29	x0	0.000000	0.000000	0.044076
	30	x0	0.000000	0.000000	0.045953
	31	x0	0.000000	0.000000	0.046945
	32	x0	0.000000	0.000000	0.049003
	33	x0	0.000000	0.000000	0.048688
	34	x0	0.000000	0.000000	0.050613
	35	x0	0.000000	0.000000	0.048918
	36	x0	0.000000	0.000000	0.048156
	37	x0	0.000000	0.000000	0.046945
	38	x0	0.000000	0.000000	0.042006
	39	x0	0.000000	0.000000	0.037709

As in MLE models, the output values of the Bayesian models represent probability density values. Note that we also have a fixed time point.

Now let’s see what happens if we reduce the distance of the stimuli:

res = model.run(auditory_position=-3, visual_position=3, causes_kind="count")
res

<NDResult 'Kording2007', modes=['auditory' 'visual' 'multi'], times=1, positions=40, positions_coordinates=1, causes=1>

The model now reports that the stimuli as originating from a common source. We can directly observe the probability of the stimuli originating from a common cause:

print("p(C=1):", res.e_["mean_p_common_cause"])
print("C:", res.causes_)

p(C=1): 0.561605322439197
C: 1

Refer to the API documentation for more details about the models available in the bayesian module.

This demonstration of the Bayesian Causal Inference model mechanics is inspired in the tutorial of the BCIT Toolbox.

Implementation of neural models

You can implement the network model developed by Cuppini et al. (2017) on by importing the corresponding module and instantiating the Cuppini2017 class:

from skneuromsi.neural import Cuppini2017

model = Cuppini2017(neurons=90, position_range=(0, 90))
model

<skneuromsi.neural._cuppini2017.Cuppini2017 at 0x708db6521180>

vars(model)

{'_neurons': 90,
 '_position_range': (0, 90),
 '_position_res': 1.0,
 '_time_range': (0, 100),
 '_time_res': 0.01,
 '_integrator': <brainpy._src.integrators.ode.explicit_rk.Euler at 0x708db6434430>,
 '_random': Generator(PCG64) at 0x708DC5052880,
 '_mode0': 'auditory',
 '_mode1': 'visual',
 'run': <function skneuromsi.neural._cuppini2017.Cuppini2017.run(*, auditory_position=None, visual_position=None, auditory_sigma=32, visual_sigma=4, auditory_intensity=28, visual_intensity=27, auditory_duration=None, auditory_onset=0, auditory_stim_n=1, visual_duration=None, visual_onset=0, visual_stim_n=1, auditory_soa=None, visual_soa=None, noise=False, noise_level=0.4, feedforward_weight=18, cross_modal_weight=1.4, causes_kind='count', causes_dim='space', causes_peak_threshold=0.15, causes_peak_distance=None)>}

You can refer to the API documentation for more details about the available parameters.

As before, let’s run the model for two conflicting stimulus locations:

res = model.run(auditory_position=30, visual_position=60)

vars(res)

{'_mname': 'Cuppini2017',
 '_mtype': 'Neural',
 '_output_mode': 'multi',
 '_nmap': {'auditory': 'auditory', 'visual': 'visual'},
 '_time_range': array([  0, 100]),
 '_position_range': array([ 0, 90]),
 '_time_res': 0.01,
 '_position_res': 1.0,
 '_run_parameters': {'auditory_position': 30,
  'visual_position': 60,
  'auditory_sigma': 32,
  'visual_sigma': 4,
  'auditory_intensity': 28,
  'visual_intensity': 27,
  'auditory_duration': None,
  'auditory_onset': 0,
  'auditory_stim_n': 1,
  'visual_duration': None,
  'visual_onset': 0,
  'visual_stim_n': 1,
  'auditory_soa': None,
  'visual_soa': None,
  'noise': False,
  'noise_level': 0.4,
  'feedforward_weight': 18,
  'cross_modal_weight': 1.4,
  'causes_kind': 'count',
  'causes_dim': 'space',
  'causes_peak_threshold': 0.15,
  'causes_peak_distance': None},
 '_extra': {'causes_kind': 'count',
  'causes_dim': 'space',
  'causes_peak_threshold': 0.15,
  'causes_peak_distance': None,
  'stim_position': [30, 60]},
 '_causes': 2,
 '_nddata': <xarray.DataArray 'values' (modes: 3, times: 10000, positions: 90,
                             positions_coordinates: 1)> Size: 22MB
 array([[[[1.19097530e-03],
          [1.31460645e-03],
          [1.44163460e-03],
          ...,
          [8.54958672e-04],
          [9.59982932e-04],
          [1.07232604e-03]],

         [[2.24271407e-03],
          [2.48409320e-03],
          [2.73379601e-03],
          ...,
          [1.59484536e-03],
          [1.79605986e-03],
          [2.01258518e-03]],

         [[3.17001717e-03],
          [3.52222241e-03],
          [3.88888969e-03],
          ...,
 ...
          ...,
          [2.31613388e-03],
          [2.33042543e-03],
          [2.33972890e-03]],

         [[2.34435053e-03],
          [2.34445919e-03],
          [2.34007850e-03],
          ...,
          [2.31613358e-03],
          [2.33042520e-03],
          [2.33972873e-03]],

         [[2.34435040e-03],
          [2.34445910e-03],
          [2.34007844e-03],
          ...,
          [2.31613328e-03],
          [2.33042497e-03],
          [2.33972856e-03]]]])
 Coordinates:
   * modes                  (modes) <U8 96B 'auditory' 'visual' 'multi'
   * times                  (times) int64 80kB 0 1 2 3 4 ... 9996 9997 9998 9999
   * positions              (positions) int64 720B 0 1 2 3 4 5 ... 85 86 87 88 89
   * positions_coordinates  (positions_coordinates) <U2 8B 'x0'}

The network models are characterized for having multiple parameters that define the network behavior. Also, these models exploit the temporal dimension since the neural activity is computed as a numerical integration at every time step:

res.get_modes()

		modes	auditory	multi	visual
times	positions	positions_coordinates
0	0	x0	0.001191	0.000025	0.000002
	1	x0	0.001315	0.000025	0.000002
	2	x0	0.001442	0.000025	0.000002
	3	x0	0.001570	0.000025	0.000002
	4	x0	0.001699	0.000025	0.000002
...	...	...	...	...	...
9999	85	x0	0.000267	0.002271	0.000015
	86	x0	0.000309	0.002296	0.000015
	87	x0	0.000356	0.002316	0.000015
	88	x0	0.000412	0.002330	0.000015
	89	x0	0.000476	0.002340	0.000016

900000 rows × 3 columns

The output values of this models represent neural activity values. Note that the model output now haves multiple time points, so we have the neural activity values at each time point of the simluation for each spatial coordinate encoded in the network.

We observe that the model detects two distinct causes from the stimuli by looking at the causes output of the model:

print("C:", res.causes_)

C: 2

Now let’s see what happens if we reduce the distance of the stimuli:

res = model.run(auditory_position=40, visual_position=50)
res

<NDResult 'Cuppini2017', modes=['auditory' 'visual' 'multi'], times=10000, positions=90, positions_coordinates=1, causes=1>

print("C:", res.causes_)

C: 1

The model outputs a common cause for stimuli close in space.

Refer to the API documentation for more details about the models available in the neural module.