Study of the cost distributions found by Bayes

>> clear
>> load info_posreal_zonapeque
>> prob=infoarchivos2prob(infoarchivos,1:6);

Parameters of the most likely value:

>> [m,ind_max]=max(prob(:));
>> [pot_max,alfa_max,beta_max,lambda_max,mu_max,tipo_max]=ind2sub(size(prob),ind_max);

Now, we will study all values with a relatively high probability, but only those with alfa, beta and pot equal to the most likely ones:

>> ind=find(prob(:)>m/100);
>> [c_pot,c_alfa,c_beta,c_lambda,c_mu,c_tipo]=ind2sub(size(prob),ind);
>> ind=ind(c_pot==pot_max & c_alfa==alfa_max & c_beta==beta_max);
>> [loglambda,logmu,tipo]=infoarchivos2parametrosdistrib(infoarchivos,ind)
>> x=0:.001:4;
>> P=parametrosdistrib2prob(x,loglambda,logmu,tipo);
>> plot(x,P)

>> semilogy(x,P)

>> loglog(x,P)

Very similar, although coming from different types:

>> tipo

 

tipo =

 

1.0000 -1.0000

2.1000 1.0000

2.1000 2.0000

2.1000 2.0000

2.2000 1.0000

2.2000 1.0000

2.2000 1.0000

2.2000 1.0000

2.2000 1.0000

Advertisement

New Figure 1

>> fig_pap_presentacion_05(1,zeros(100))

New version:

New Figure 2

>> fig_pap_elegans_15(1,zeros(100))

I have added the text “Simulation” and “Experiment” Do you like it?

I change it after Comment 1:

>> fig_pap_elegans_15(1,zeros(100))

I change it to remove “predicted”:

New Figure 3

>> fig_pap_bayes_09(1,zeros(100))

Only the axislabels of box A and the legend of box B have changed.

Any idea to improve the axislabels?

The reviewer asked to write “that $\alpha=\beta=1/29.3$ corresponds to the ideal wiring economy”. But I thought that ‘ideal wiring economy values’ was too long. Is it correct this way?

New Figure 4 after revision

>> fig_pap_coli_10(1,zeros(100))

Main questions:

– Are the labels OK? In particular, this style with the description, comma, the symbol. And the description of Delta.
– Are the words ATP, BIOMASS where they should be? In boxes A, B they might be inside the box.
– Colorbar: It is modified to show dark red in the top. Enough?

Saving information about known networks

>> chemotaxis1.info=’Red de chemotaxis que aparece en el primer post sobre dinámica de elegans’;
>> chemotaxis1.buenas=buenas;
>> chemotaxis1.buenas=[11:14 19:20 40:41 54:57 77:78 111:113 168:169 178:179];
>> chemotaxis1.sensoriales=[40 41];
>> chemotaxis1.motoras=54:57;
>> save info_redes chemotaxis1

 

Comparison of different dynamic models for C elegans

I use something like a ROC curve, that measures the proportion of the neurons belonging to the subnetwork and the proportion of neurons not belonging to the subnetwork that the method finds, depending on the threshold:

I use the chemotaxis network I used before.

Comparison between the model with all conductances equal (as in the older post) and the model with all connections as if gap junctions, and no membrane conductance:

>> clear
>> load(‘c:\hipertec\ablationelegans\Datos\datos_celegans_paradinamica.mat’)
>> V0=zeros(279,1);
>> V0([40 41])=1;
>> M{1}=conectividad2matrizsistema(todas.A_chem,todas.A_ej,todas.GABA,1,1,1);
>> resultados{1}=matrizsistema2experimentoablacion(M{1},V0,500,.005,[-1 1],54:57);
>> M{2}=conectividad2matrizsistema(zeros(279),todas.A_chem+todas.A_ej,todas.GABA,1,1,0);
>> resultados{2}=matrizsistema2experimentoablacion(M{2},V0,500,.005,[-1 1],54:57);
>> buenas=[11:14 19 20 77 78 111:113 168 169 178 179]; % Neurons belonging to the subnetwork
>> resultados2pseudoroc(resultados{1},buenas,[40 41 54:57],’b’);
>> hold on
>> resultados2pseudoroc(resultados{2},buenas,[40 41 54:57],’r’);

 

Now I compare the old model (with all conductances equal) with different time scales:

>> resultados{1,2}=matrizsistema2experimentoablacion(M{1},V0,250,.005,[-1 1],54:57);
>> resultados{1,3}=matrizsistema2experimentoablacion(M{1},V0,100,.005,[-1 1],54:57);
>> clf
>> resultados2pseudoroc(resultados{1,1},buenas,[40 41 54:57],’b’);
>> hold on
>> resultados2pseudoroc(resultados{1,2},buenas,[40 41 54:57],’r’);
>> resultados2pseudoroc(resultados{1,3},buenas,[40 41 54:57],’k’);

Perfect overlapping, so the time scale does not make any effect in this range.

Test of the pseudoROC:

>> resultados=rand(279,1);
>> clf
>> resultados2pseudoroc(resultados,buenas,[40 41 54:57],’k’);

OK.

Now, I compare the model with equal conductances with a model tuned so that it does not saturate:

>> M{1}=conectividad2matrizsistema(todas.A_chem,todas.A_ej,todas.GABA,1,1,1);
>> resultados{1}=matrizsistema2experimentoablacion(M{1},V0,500,.005,[-1 1],54:57);
>> resultados{1}=matrizsistema2experimentoablacion(M{1},V0,500,.005,[-1 1],54:57);
>> M{2}=conectividad2matrizsistema(todas.A_chem,todas.A_ej,todas.GABA,.05,1,1.5);
>> V=matrizsistema2V_num(M{2},V0,10000,.005,[-1 1]);
>> plot(V(54:57,:)’)

>> resultados{2}=matrizsistema2experimentoablacion(M{2},V0,10000,.005,[-1 1],54:57);
>> resultados2pseudoroc(resultados{2},buenas,[40 41 54:57],’r’);
>> hold on
>> resultados2pseudoroc(resultados{1},buenas,[40 41 54:57],’b’);

The model with equal conductances works better.

First results for network identification

I just run the network with the linear model with saturation, and all conductances equal to 1. Then, I ablate neurons one by one, and see the change in the target neurons. I am looking for the following network:

>> M=conectividad2matrizsistema(todas.A_chem,todas.A_ej,todas.GABA,1,1,1);

Initial conditions, the sensory neurons activated:

>> find(strcmpi(todas.nombres,’aser’))
ans =
41
>> find(strcmpi(todas.nombres,’asel’))
ans =
40
>> V0=zeros(279,1);
>> V0([40 41])=1;

Target neurons:

>> find(strcmpi(todas.nombres,’avar’))
ans =
55
>> find(strcmpi(todas.nombres,’aval’))
ans =
54
>> find(strcmpi(todas.nombres,’avbl’))
ans =
56
>> find(strcmpi(todas.nombres,’avbr’))
ans =
57
>> resultados=matrizsistema2experimentoablacion(M,V0,500,.005,[-1 1],54:57);
>> plot(resultados)

The highest peaks are for the target neurons. Zoom:

The high peaks around 40 are the target neurons, but the other high peak are other neurons:

>> todas.nombres(find(mean(resultados,2)>4*10^-3 & mean(resultados,2)<10^-2))
ans =
‘AIBL’
‘AIBR’
‘ASEL’
‘ASER’

AIB, that plays a central role in the network. As we lower the threshold, we find more and more neurons in the network, and also some which are not:

>> todas.nombres(find(mean(resultados,2)>1*10^-3 & mean(resultados,2)<10^-2))
ans =
‘AIBL’
‘AIBR’
‘AIYL’
‘AIYR’
‘AIZL’
‘AIZR’
‘ASEL’
‘ASER’
‘AVEL’
‘AVER’
‘PVCL’
‘PVCR’
‘RIBL’
‘RIML’
‘RIMR’

We catch 7 out of 15 neurons in the network (counting laterally doubled neurons, excluding sensory and target neurons). 6 neurons are not in the network.

>> todas.nombres(find(mean(resultados,2)>.5*10^-3 & mean(resultados,2)<10^-2))
ans =
‘AIAL’
‘AIAR’
‘AIBL’
‘AIBR’
‘AIYL’
‘AIYR’
‘AIZL’
‘AIZR’
‘ASEL’
‘ASER’
‘AVDR’
‘AVEL’
‘AVER’
‘DVA’
‘FLPL’
‘PVCL’
‘PVCR’
‘RIBL’
‘RIBR’
‘RIML’
‘RIMR’
‘SAAVL’
‘VA08’

We catch 11 out of 15 neurons. 10 neurons are not in the network.

>> todas.nombres(find(mean(resultados,2)>.4*10^-3 & mean(resultados,2)<10^-2))
ans =
‘AIAL’
‘AIAR’
‘AIBL’
‘AIBR’
‘AIYL’
‘AIYR’
‘AIZL’
‘AIZR’
‘AS07’
‘ASEL’
‘ASER’
‘AVDL’
‘AVDR’
‘AVEL’
‘AVER’
‘AWCL’
‘DVA’
‘DVC’
‘FLPL’
‘PVCL’
‘PVCR’
‘RIBL’
‘RIBR’
‘RIFR’
‘RIML’
‘RIMR’
‘SAAVL’
‘VA08’
‘VB09’

We catch 13 out of 15 neurons. 14 neurons are not in the network.

 

Carpenter plots of the bayesian estimator

>> t_reaccion=bayestim_medias_nosave(10^5,10,1000,2,0,1,10^-3,.995);
>> hist(t_reaccion,0:1000)

>> t_reaccion2histog_carpenter(t_reaccion,20,’r.-‘);

Mistake creating the system matrix

In previous posts, I made a mistake: I used the transposed of the chemical connectivity matrix. Fortunately, it seems there is not a fundamental difference:

Transposed,

>> M=conectividad2matrizsistema(todas.A_chem’,todas.A_ej,todas.GABA,1,1,30);
>> solucion=matrizsistema2solucion(M,V0);
>> plot(solucion.autovalores,’.’)

Not transposed:

>> M=conectividad2matrizsistema(todas.A_chem,todas.A_ej,todas.GABA,1,1,30);
solucion=matrizsistema2solucion(M,V0);
plot(solucion.autovalores,’.’)