Demo¶

This tutorial is to explain how we can go from simulation to estimation with nm-bridge and vice-versa.

First we need to load the libraries (nm-bridge, the main one + igraph to plot graphs) and a source code for making actually the plots

In [1]:

library(nmbridge)
library(igraph)
source('plotting.R')

library(nmbridge) library(igraph) source('plotting.R')

Warning message:
“package ‘nmbridge’ was built under R version 4.3.0”

Attaching package: ‘igraph’


The following objects are masked from ‘package:stats’:

    decompose, spectrum


The following object is masked from ‘package:base’:

    union

Simple Parameters for the simulations¶

nm-bridge allows to simulate complex patterns : switch in behaviors, multiple trials, with potentially different lengths. The simple parameters consists in

• Ntrial, the number of trials
• Ncomp, the number of different behaviors that you want
• Nneur, the number of neurons involved
• delay, the time before actually collecting the data : Hawkes processes need to warm up before reaching stationarity. This warm up period is at the beginning of each trial, it is not in between behaviors inside a simulation.

In [2]:

delay = 0 
Nneur = 3 
Ntrial= 4 
Ncomp = 2

delay = 0 Nneur = 3 Ntrial= 4 Ncomp = 2

Behaviors¶

Each behavior $m$ corresponds to a set of spontaneous parts (one per neuron $i$) and a set of interaction functions from $j$ to $i$ for every other neuron $j$ (including $i$ itself, this is the self-interaction part). So the intensity of the spike train $N^i$ of neuron $i$ during behavior $m$ (that is time $t$ is in a period where behavior is $m$) is given by $$\lambda^i_t=\nu^i_m +\int_{-\infty}^t h^{j\to i}_m(t-s) dN^j_s$$ Note that when at the boundary between two behaviors the past in the integral is in made of times $s$, that might belong to another behavior.

the spontaneous matrix

We put all the $\nu^i_m$ in a matrix, where the column represents the behavior $m$, and the line the neuron $i$.

In [3]:

spontComp = matrix(c(20,10,3,30,40,0), nrow=3) 

spontComp

spontComp = matrix(c(20,10,3,30,40,0), nrow=3) spontComp

A matrix: 3 × 2 of type dbl

20	30
10	40
3	0

In this example for instance, Neuron 3 has a null spontaneous part in the second behavior.

the array of interaction functions

We put all the $h^{j\to i}_m(.)$ in an array of size 3 of cells.

• First index is $i$ the target neuron (to which ?),
• second index is $j$ the targetting neuron (from which ?)
• third index is $m$ the behavior.

Once the array is created, at each position a matrix with two columns gives the piecewise constant representation of the function $h^{j\to i}_m(.)$

In [4]:

interacComp= array(c(list()),c(Nneur,Nneur,Ncomp)) # 1st index to which, 2nd from which, 3rd behavior

# Behavior 1 
interacComp[1,1,1][[1]]= matrix(c(0, 0.02, 30, 0), nrow = 2)
interacComp[2,1,1][[1]]= matrix(c(0, 0.01, 0.02, 30, 0, 0), nrow = 3)
interacComp[3,1,1][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2)

interacComp[1,2,1][[1]]= matrix(c(0, 0.01, 0.02, 0, 30, 0), nrow = 3)
interacComp[2,2,1][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2)
interacComp[3,2,1][[1]]= matrix(c(0, 0.02, 10, 0), nrow = 2)

interacComp[1,3,1][[1]]= matrix(c(0, 0.01, 0.02, 0, 30, 0), nrow = 3)
interacComp[2,3,1][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2)
interacComp[3,3,1][[1]]= matrix(c(0, 0.02, 10, 0), nrow = 2)

# Behavior 2

interacComp[1,1,2][[1]]= matrix(c(0, 0.01, 0.02, 0, 30, 0), nrow = 3)
interacComp[2,1,2][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2)
interacComp[3,1,2][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2)

interacComp[1,2,2][[1]]= matrix(c(0, 0.02, 30, 0), nrow = 2)
interacComp[2,2,2][[1]]= matrix(c(0, 0.01, 0.02, 30, 0, 0), nrow = 3)
interacComp[3,2,2][[1]]= matrix(c(0, 0.02, 1, 0), nrow = 2)

interacComp[1,3,2][[1]]= matrix(c(0, 0.02, 30, 0), nrow = 2)
interacComp[2,3,2][[1]]= matrix(c(0, 0.01, 0.02, 30, 0, 0), nrow = 3)
interacComp[3,3,2][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2)

interacComp[1,2,1][[1]]

interacComp= array(c(list()),c(Nneur,Nneur,Ncomp)) # 1st index to which, 2nd from which, 3rd behavior # Behavior 1 interacComp[1,1,1][[1]]= matrix(c(0, 0.02, 30, 0), nrow = 2) interacComp[2,1,1][[1]]= matrix(c(0, 0.01, 0.02, 30, 0, 0), nrow = 3) interacComp[3,1,1][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2) interacComp[1,2,1][[1]]= matrix(c(0, 0.01, 0.02, 0, 30, 0), nrow = 3) interacComp[2,2,1][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2) interacComp[3,2,1][[1]]= matrix(c(0, 0.02, 10, 0), nrow = 2) interacComp[1,3,1][[1]]= matrix(c(0, 0.01, 0.02, 0, 30, 0), nrow = 3) interacComp[2,3,1][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2) interacComp[3,3,1][[1]]= matrix(c(0, 0.02, 10, 0), nrow = 2) # Behavior 2 interacComp[1,1,2][[1]]= matrix(c(0, 0.01, 0.02, 0, 30, 0), nrow = 3) interacComp[2,1,2][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2) interacComp[3,1,2][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2) interacComp[1,2,2][[1]]= matrix(c(0, 0.02, 30, 0), nrow = 2) interacComp[2,2,2][[1]]= matrix(c(0, 0.01, 0.02, 30, 0, 0), nrow = 3) interacComp[3,2,2][[1]]= matrix(c(0, 0.02, 1, 0), nrow = 2) interacComp[1,3,2][[1]]= matrix(c(0, 0.02, 30, 0), nrow = 2) interacComp[2,3,2][[1]]= matrix(c(0, 0.01, 0.02, 30, 0, 0), nrow = 3) interacComp[3,3,2][[1]]= matrix(c(0, 0.02, 0, 0), nrow = 2) interacComp[1,2,1][[1]]

A matrix: 3 × 2 of type dbl

0.00	0
0.01	30
0.02	0

In the above example, $h^{2\to1}_{1}$ is equal to 0 between 0 and 0.01, then 30 between 0.01 and 0.02 and is null afterwards. One always needs to have a 0 in the the top left and bottom right of these functions.

Matrix of composition of each trial in terms of behavior¶

We need to give the succession of behaviors to the simulator. To do that we need a matrix with 4 columns

• First column : the index of the trial
• Second column : the beginning of a period (box) with a certain behavior
• Third column : the end of this period (box)
• Fourth column : the index of the behavior during the period

In [5]:

BoxComp = t(matrix(c(1, 0, 50, 1, 1, 50, 100, 2, 2, 0, 15, 2, 2, 15, 35,1,3,0,10,1,4,0,40,2), nrow=4) )

BoxComp

BoxComp = t(matrix(c(1, 0, 50, 1, 1, 50, 100, 2, 2, 0, 15, 2, 2, 15, 35,1,3,0,10,1,4,0,40,2), nrow=4) ) BoxComp

A matrix: 6 × 4 of type dbl

1	0	50	1
1	50	100	2
2	0	15	2
2	15	35	1
3	0	10	1
4	0	40	2

Here for instance, the first trial is very long (100 s) with 2 behaviors, whereas the third trial (10s) is very short with only the first behavior.

Be careful that the boxes should always form a partition of the interval 0 to the length of the trials. If there are holes or if they overlap, the simulation will not work.

If delay is positive, before the beginning of each trial, the simulator will run the simulation with the first encountered behavior of the trial on $[-delay, 0]$, these data are not saved but just help to give a non empty past when starting in 0. Before -delay, the past is empty.

If delay is null, the simulation starts with an empty past before 0.

Simulation¶

There are two functions to simulate depending on the form of output that you want.

The DataNeur shape

In this shape, the data are put inside a matrix, with as many lines as $Nneur*Ntrial$. Each line corresponds to one neuron in one trial (first all trials of neuron 1, then all trials of neuron 2 etc). The first element of the line is the total number of spikes of this neuron in this trial, then we have the times of the spikes in increasing order. These times are always positive, because BoxComp is always proposing non negative time for the boxes.

The number of columns is the maximal number of spikes per neuron and trials +1.

For the neurons/trials who do not have this maximal number of spikes, the line is filled with 0 after the last spike.

The parameter seeds is to fix the random seed in the simulation and help reproducibility of the results.

Note that the function (as the other one) do not need Nneur and Ntrial, because these informations are included in the shape of spontComp for instance.

In [6]:

dataDn = HawkesmultiDn(Ntrial, delay, BoxComp, spontComp, interacComp, seeds=5)

dataDn[,0:15]

dataDn = HawkesmultiDn(Ntrial, delay, BoxComp, spontComp, interacComp, seeds=5) dataDn[,0:15]

A matrix: 12 × 15 of type dbl

8902	0.00416614	0.01147903	0.01546993	0.01859347	0.01887805	0.01915638	0.03083862	0.03324425	0.04597204	0.04801024	0.04808729	0.04842853	0.04924836	0.05481737
3183	0.01351614	0.02853352	0.03128008	0.04215050	0.10097411	0.12177982	0.12298958	0.12534916	0.13277432	0.13740138	0.18840377	0.19291055	0.20410897	0.22066814
884	0.05168480	0.05690939	0.05745177	0.07497112	0.09650471	0.09940522	0.11437276	0.12706292	0.24042556	0.24494447	0.24673938	0.29151798	0.38217403	0.44691457
3631	0.02115912	0.04511300	0.13108487	0.13308240	0.16181567	0.16402041	0.17044417	0.18142204	0.18622657	0.19329798	0.21067264	0.21495744	0.22063040	0.22590757
4792	0.03365308	0.04264822	0.05651689	0.06055358	0.08022693	0.09317943	0.10187255	0.10416393	0.11379479	0.12003039	0.12193579	0.12650694	0.13386749	0.14615059
1561	0.00878259	0.01116875	0.03079931	0.11543262	0.12075950	0.12983041	0.17268687	0.18526586	0.20158794	0.23795363	0.25650694	0.26305924	0.29164201	0.29437908
338	0.08008036	0.10894502	0.11650925	0.15254469	0.33041388	0.34903977	0.45153624	0.45441338	0.45863711	0.47324692	0.47556022	0.48571867	0.48765381	0.48873784
2215	0.01207430	0.17058410	0.20521414	0.22265121	0.22555818	0.22625960	0.22649720	0.23513593	0.27894897	0.28313858	0.29223468	0.29458475	0.30049764	0.31063532
669	0.03485429	0.05062657	0.05162470	0.05452395	0.06646709	0.15295788	0.23798365	0.24132308	0.25810334	0.37630004	0.43818581	0.63479173	0.73028808	0.74403311
283	0.13842273	0.55544235	3.08905757	3.44240568	3.44969486	4.05025883	4.94914798	4.95673032	5.29004787	6.54139090	9.10834082	9.69382020	9.74002280	11.00983348
125	0.21147630	0.34220476	0.50914456	0.52255259	0.56181063	0.56795014	0.56858097	0.58296928	0.59041783	0.60255365	0.76838272	0.77399834	0.77538418	1.04552682
44	4.40341703	4.53645578	5.99730298	6.45519158	7.82258614	8.15700570	8.16336071	10.84365874	10.91747656	13.47800051	13.81040656	14.15211144	14.34124097	15.03987445

Note that we have always a very small number of spikes in trial 3 (whatever the neuron) and that the most active neuron whatever the trial seems to be neuron 1.

We can also visualise the spikes.

In [7]:

spikes.plot(dataDn, n=Nneur, col=c('red','green','blue'))

spikes.plot(dataDn, n=Nneur, col=c('red','green','blue'))

Be careful, trials are increasing vertically, so the one at the bottom is the first trial. One can also restrict the time window

In [8]:

spikes.plot(dataDn, n=Nneur,xrange=c(0,10), col=c('red','green','blue'))

spikes.plot(dataDn, n=Nneur,xrange=c(0,10), col=c('red','green','blue'))

Remember that in Behavior 2, neuron 3 should be less active because it has a null spontaneous activity (and (also the interaction functions are small). You can see that in trials 2 and 4 (from the bottom)

The Data list shape

Here it produces a list of size Ntrial. Each element of the list corresponds to a trial and it is a two column matrix, where

• the first column contains all the spikes of all neurons in this trial in increasing order
• the second column contains the index (minus 1) of the neuron that emitted the spike on the same line. Be careful neuron 1, correspond to index 0 in this representation (cf Python or C++ rule)

In [9]:

dataDs= HawkesmultiDs(Ntrial, delay, BoxComp, spontComp, interacComp, seeds=5)

length(dataDs)

dataDs= HawkesmultiDs(Ntrial, delay, BoxComp, spontComp, interacComp, seeds=5) length(dataDs)

4

In [10]:

dataDs[[1]][1:20,]

dataDs[[1]][1:20,]

A matrix: 20 × 2 of type dbl

0.00416614	0
0.01147903	0
0.01546993	0
0.01859347	0
0.01887805	0
0.01915638	0
0.03083862	0
0.03324425	0
0.03365308	1
0.03485429	2
0.04264822	1
0.04597204	0
0.04801024	0
0.04808729	0
0.04842853	0
0.04924836	0
0.05062657	2
0.05162470	2
0.05452395	2
0.05481737	0

We see typically that neuron 1 is spiking much more than the other ones.

Reconstruction¶

Now we want to use the data we just simulated to reconstruct the connectivity graph.

How to apply reconstruction

• k is the number of bins/pieces to reconstruct the interaction functions
• delta is the length of these bins
• Ntrial is still the number of trials (necessary to disambiguate from a DataNeur shape between Nneur and Ntrial)
• scale : the spikes are going to be rounded at $10^{-scale}$ to gain in rapidity of computation, by default it is 4
• gamma : vector of possible tuning parameters for the Bernstein lasso (3.5 is a classic choice)
• BoxComp : you need to tell the estimator where are the behaviors, for the moment we use the one above, but you could change. Here you do not need the boxes to cover the whole space.
• perBox : an indicator if TRUE, the estimation is made Box per Box if not, all the boxes with the same behavior are aggregated to get an estimator per behavior.
• colnums : an indicator that is there to give the index of the label columns and spikes in the case of a Data list shape (useful when recorded data with various formats)

In [11]:

k=2
delta=0.01
Ntrial=4
gamma=3.5

k=2 delta=0.01 Ntrial=4 gamma=3.5

In [12]:

estimDn=Lasso(dataDn,BoxComp,k,delta,perBox=FALSE,gamma,Ntrial) 
# or Lasso(dataDn,BoxComp,k,delta,perBox=FALSE,gamma,Ntrial,scale=4), this is the same

estimDn=Lasso(dataDn,BoxComp,k,delta,perBox=FALSE,gamma,Ntrial) # or Lasso(dataDn,BoxComp,k,delta,perBox=FALSE,gamma,Ntrial,scale=4), this is the same

In estimDn, you have

• estimDn$BL the Bernstein Lasso estimator. This is a list: each element = a matrix ($(Nneur*K+1)*Nneur$) corresponding to one behavior. The matrix can be read as followed : coefficients for a given neuron $i$ are in the corresponding column $i$. Then first line = spontaneous parts, next lines are coefficients of $h^{1\to i}_m$, until $h^{Nneur\to i}_m$
• estimDn$BOL the Ordinary least square estimator on the support given by the Lasso (same convention as above). This estimator should be less biased than the first one.
• estimDn$b, estimDn$g, estimDn$mu2 and estimDn$muinfty are all the matrices necessary to compute the Lasso criterion
• estimDn$msg is here to tell if there has been any trouble during the computation

In [13]:

estimDn$BOL

estimDn$BOL

1. A matrix: 7 × 3 of type dbl

   21.32501	8.93155	2.689272
   27.53378	31.60402	0.000000
   30.92597	0.00000	0.000000
   0.00000	0.00000	8.626549
   29.19262	0.00000	10.664508
   0.00000	0.00000	11.521035
   29.72309	0.00000	11.600578

2. A matrix: 7 × 3 of type dbl

   30.75399	39.16664	-0.1062111
   0.00000	0.00000	0.0000000
   29.52956	0.00000	0.0000000
   28.64010	28.32208	1.1289429
   30.66712	2.95349	1.0130973
   0.00000	0.00000	0.0000000
   0.00000	0.00000	0.0000000

In [14]:

estimDn$msg

estimDn$msg

'Successful computation for all the neurons (Nneur = 3), for any gamma in [3.5] and all the comportments (Ncomp = 2)'

Same thing is possible from a Data List shape. Just note that with colnums you have to tell the index of the column with labels and next the index of the column with spike times. The result has exactly the same shape as above.

In [15]:

colnums=c(2,1)
estimDs=Lasso(dataDs,BoxComp,k,delta,perBox=FALSE,gamma,Ntrial,colnums)

estimDs$BOL
estimDs$msg

colnums=c(2,1) estimDs=Lasso(dataDs,BoxComp,k,delta,perBox=FALSE,gamma,Ntrial,colnums) estimDs$BOL estimDs$msg

1. A matrix: 7 × 3 of type dbl

   21.20970	9.061071	2.621332
   27.95819	31.464341	0.000000
   30.66235	0.000000	0.000000
   0.00000	0.000000	8.600724
   29.32506	0.000000	10.638605
   0.00000	0.000000	11.480248
   29.41392	0.000000	11.725928

2. A matrix: 7 × 3 of type dbl

   30.75399	39.16664	-0.1062111
   0.00000	0.00000	0.0000000
   29.52956	0.00000	0.0000000
   28.64010	28.32208	1.1289429
   30.66712	2.95349	1.0130973
   0.00000	0.00000	0.0000000
   0.00000	0.00000	0.0000000

'Successful computation for all the neurons (Nneur = 3), for any gamma in [3.5] and all the comportments (Ncomp = 2)'

How to represent the reconstruction

First thing is to convert the coefficients into the same format as the matrix of spontaneous activity (S) and the array of interaction functions(I). Everything is wrapped up in a list with S for the spontaneous activities and with I for the interaction functions

In [16]:

BOL=Lasso_func(estimDn$BOL,k,delta)
BL=Lasso_func(estimDn$BL,k,delta)
True=list(S=spontComp,I=interacComp)

BOL=Lasso_func(estimDn$BOL,k,delta) BL=Lasso_func(estimDn$BL,k,delta) True=list(S=spontComp,I=interacComp)

In [17]:

BOL$S

BOL$S

A matrix: 3 × 2 of type dbl

21.325011	30.7539889
8.931550	39.1666449
2.689272	-0.1062111

In [18]:

BOL$I[1,2,1][[1]]

BOL$I[1,2,1][[1]]

A matrix: 3 × 2 of type dbl

0.00	0.00000
0.01	29.19262
0.02	0.00000

Next we can plot all the results to compare them.

• names contains the names of the behaviors
• neurnames if given, gives the names of the neurons
• give as many colors as the size of the list

In [19]:

listIS=list(True,BL,BOL)

Plot_Hawkes(listIS,c('black','green','red'),names=c('comp1','comp2'))

listIS=list(True,BL,BOL) Plot_Hawkes(listIS,c('black','green','red'),names=c('comp1','comp2'))

True values are represented in black, while BL and BOL estimators are represented in green and red respectively.

Most of the time the BOL estimator is closer than the BL.

You can also represent graphs.

• takes in entries a list made of S and I
• names and neurnames as above
• you can either not give a layout and one will be automatically given, or you force it (easier to compare graphs) (layout =position of the nodes, first column the 'x' second column the 'y')
• vecpres gives the number of lines and columns for the plots of the different behaviors

The thickness of the arrows is proportional to the L1 norm of $h^{j\to i}_m$.

In [20]:

Plot_graphs_Hawkes(True,names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2))
Plot_graphs_Hawkes(BOL,names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2))

Plot_graphs_Hawkes(True,names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2)) Plot_graphs_Hawkes(BOL,names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2))

In [21]:

layout=matrix(c(0,0.5,1,0,0.5,0),ncol=2)

layout

layout=matrix(c(0,0.5,1,0,0.5,0),ncol=2) layout

A matrix: 3 × 2 of type dbl

0.0	0.0
0.5	0.5
1.0	0.0

In [22]:

Plot_graphs_Hawkes(True,names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2),layout=layout)
Plot_graphs_Hawkes(BOL,names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2),layout=layout)

Plot_graphs_Hawkes(True,names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2),layout=layout) Plot_graphs_Hawkes(BOL,names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2),layout=layout)

Remark that in behavior 2, 3 is not very active, that is why the estimations around node 3 are not as good as they should be.

Start with real data¶

In [23]:

download.file('https://plmbox.math.cnrs.fr/f/b2b22315305e45ad9f85/?dl=1','SDM113_session9_numbers_ms.csv')

download.file('https://plmbox.math.cnrs.fr/f/b2b22315305e45ad9f85/?dl=1','SDM113_session9_numbers_ms.csv')

In practice it is more likely that we will start with real data. You can download a small subpart of the striatum/free learning task experiment, here. Data are in milliseconds so for interpretation it is easier if we convert them.

In [24]:

ds = as.matrix(read.table("SDM113_session9_numbers_ms.csv",sep=',',header=T))
ds[,2]=ds[,2]/1000 # to go to seconds

ds = as.matrix(read.table("SDM113_session9_numbers_ms.csv",sep=',',header=T)) ds[,2]=ds[,2]/1000 # to go to seconds

In [25]:

ds[1:10,]

ds[1:10,]

A matrix: 10 × 2 of type dbl

process	time
1	0.0339
1	0.0364
1	0.0437
1	0.0821
1	0.1121
1	0.1272
1	0.1313
1	0.1367
1	0.1462
1	0.1532

Here ds is in the form of a Data shape list, except that the first column is the index, the second the spike times

The annotation for two different paths is here.

In [26]:

download.file('https://plmbox.math.cnrs.fr/f/1fbb48b46d3e4738b165/?dl=1','SDM113_session9_comportment_numbers.txt')

download.file('https://plmbox.math.cnrs.fr/f/1fbb48b46d3e4738b165/?dl=1','SDM113_session9_comportment_numbers.txt')

In [27]:

boxes = as.matrix(read.table("SDM113_session9_comportment_numbers.txt",sep=',', header=T))
boxes[1:10,]

boxes = as.matrix(read.table("SDM113_session9_comportment_numbers.txt",sep=',', header=T)) boxes[1:10,]

A matrix: 10 × 4 of type dbl

trial	begin	end	behavior
1	67.063	74.836	2
1	120.867	130.271	1
1	148.517	161.163	1
1	168.464	173.345	2
1	195.370	201.256	2
1	205.788	210.444	1
1	220.123	223.364	2
1	229.886	233.742	1
1	237.556	241.065	2
1	245.204	248.404	1

It is the same format as a BoxComp, the matrix of composition of the trials, except that here we have only one trial, but a lot of different boxes. They do not cover the whole recording but it does not matter for the reconstruction part. We can apply the Lasso estimator.

In [28]:

k = 2
delta = 0.01
scale = 4 
gamma = 3.5
colnums = c(1,2)
perBox=FALSE

estimComp = Lasso(ds, boxes, k, delta, perBox, gamma, NULL, colnums, scale=scale)

estimComp$BOL

estimComp$msg

k = 2 delta = 0.01 scale = 4 gamma = 3.5 colnums = c(1,2) perBox=FALSE estimComp = Lasso(ds, boxes, k, delta, perBox, gamma, NULL, colnums, scale=scale) estimComp$BOL estimComp$msg

1. A matrix: 29 × 14 of type dbl

   15.50167	2.039568	3.975352	0.529002	0	0	0	0.4951459	1.409261	0	0.2369929	0	0.787155	0.253921
   19.74992	2.400887	1.828686	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   15.60387	1.318734	2.590985	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   16.93313	12.819184	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	10.843276	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   33.22555	0.000000	9.287811	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	20.391228	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000
   0.00000	0.000000	0.000000	0.000000	0	0	0	0.0000000	0.000000	0	0.0000000	0	0.000000	0.000000

2. A matrix: 29 × 14 of type dbl

   18.77994	2.74837	1.917366	1.414109	0	0	0	0.4200737	0.6613041	0	0	0	0.598917	0
   20.62139	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   13.04860	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	15.39018	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   29.09400	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0
   0.00000	0.00000	0.000000	0.000000	0	0	0	0.0000000	0.0000000	0	0	0	0.000000	0

'(numbering of comportments and neurons starts at 0!), g is not invertible for comportement of index 1, '

It remains to visualise it.

In [29]:

IS_data = Lasso_func(estimComp$BOL, k, delta)

IS_data = Lasso_func(estimComp$BOL, k, delta)

In [30]:

Plot_graphs_Hawkes(IS_data, names=c('comp1','comp2'), neurnames=seq(1,14), vecpres=c(1,2))

Plot_graphs_Hawkes(IS_data, names=c('comp1','comp2'), neurnames=seq(1,14), vecpres=c(1,2))

Let us now try to simulate with the parameters, and to see something, let us focus on neurons 1, 2 and 3.

In [31]:

IS_res=list(S=IS_data$S[1:3,],I=IS_data$I[1:3,1:3,])

IS_res=list(S=IS_data$S[1:3,],I=IS_data$I[1:3,1:3,])

In [32]:

Plot_graphs_Hawkes(IS_res, names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2))

Plot_graphs_Hawkes(IS_res, names=c('comp1','comp2'),neurnames=1:3, vecpres=c(1,2))

We now have to define BoxComp. For good visualisation, I make 25 trials in comp1 and 25 trials in comp2

In [33]:

Ntrial=50
BoxComp=matrix(0,ncol=4,nrow=50)
for(n in 1:25)
{
    BoxComp[n,]=c(n,0,6,1)
    BoxComp[n+25,]=c(n+25,0,6,2)
}
delay=0
spontComp=IS_res$S
interacComp=IS_res$I
BoxComp

Ntrial=50 BoxComp=matrix(0,ncol=4,nrow=50) for(n in 1:25) { BoxComp[n,]=c(n,0,6,1) BoxComp[n+25,]=c(n+25,0,6,2) } delay=0 spontComp=IS_res$S interacComp=IS_res$I BoxComp

A matrix: 50 × 4 of type dbl

1	0	6	1
2	0	6	1
3	0	6	1
4	0	6	1
5	0	6	1
6	0	6	1
7	0	6	1
8	0	6	1
9	0	6	1
10	0	6	1
11	0	6	1
12	0	6	1
13	0	6	1
14	0	6	1
15	0	6	1
16	0	6	1
17	0	6	1
18	0	6	1
19	0	6	1
20	0	6	1
21	0	6	1
22	0	6	1
23	0	6	1
24	0	6	1
25	0	6	1
26	0	6	2
27	0	6	2
28	0	6	2
29	0	6	2
30	0	6	2
31	0	6	2
32	0	6	2
33	0	6	2
34	0	6	2
35	0	6	2
36	0	6	2
37	0	6	2
38	0	6	2
39	0	6	2
40	0	6	2
41	0	6	2
42	0	6	2
43	0	6	2
44	0	6	2
45	0	6	2
46	0	6	2
47	0	6	2
48	0	6	2
49	0	6	2
50	0	6	2

In [34]:

dataDn = HawkesmultiDn(Ntrial,delay, BoxComp, spontComp, interacComp,seeds=0)

dataDn = HawkesmultiDn(Ntrial,delay, BoxComp, spontComp, interacComp,seeds=0)

In [35]:

spikes.plot(dataDn, n=3,col=c('red','green','blue'))

spikes.plot(dataDn, n=3,col=c('red','green','blue'))