diff --git a/VBA_BMA.m b/VBA_BMA.m index 5e61fca0..36b8867d 100644 --- a/VBA_BMA.m +++ b/VBA_BMA.m @@ -42,6 +42,7 @@ % hidden states T = size(p0{1}.muX,2); +try for t=1:T mus = cell(K,1); Qs = cell(K,1); @@ -51,36 +52,41 @@ end [p_BMA.muX(:,t),p_BMA.SigmaX.current{t}] = get2moments(mus,Qs,ps); end +end % data precision mus = cell(K,1); Qs = cell(K,1); -for k=1:K - mus{k} = p0{k}.a_sigma/p0{k}.b_sigma; - Qs{k} = p0{k}.a_sigma/p0{k}.b_sigma^2; +for iS=1:numel(p0{1}.a_sigma) % loop over sources + for k=1:K + mus{k} = p0{k}.a_sigma(iS)/p0{k}.b_sigma(iS); + Qs{k} = p0{k}.a_sigma(iS)/p0{k}.b_sigma(iS)^2; + end + [m,v] = get2moments(mus,Qs,ps); + p_BMA.b_sigma(iS) = m/v; + p_BMA.a_sigma(iS) = m*p_BMA.b_sigma(iS); end -[m,v] = get2moments(mus,Qs,ps); -p_BMA.b_sigma = m/v; -p_BMA.a_sigma = m*p_BMA.b_sigma; % hidden state precision mus = cell(K,1); Qs = cell(K,1); id = zeros(K,1); + for k=1:K mus{k} = p0{k}.a_alpha/p0{k}.b_alpha; Qs{k} = p0{k}.a_alpha/p0{k}.b_alpha^2; - if isinf(p0{k}.a_alpha) && p0{k}.b_alpha==0 + if (isempty(p0{k}.a_alpha) && isempty(p0{k}.b_alpha)) || (isinf(p0{k}.a_alpha) && p0{k}.b_alpha==0) id(k) = 1; end end + if isequal(sum(id),K) % all deterministic systems p_BMA.b_alpha = Inf; p_BMA.a_alpha = 0; elseif isequal(sum(id),0) % all stochastic systems [m,v] = get2moments(mus,Qs,ps); p_BMA.b_alpha = m/v; - p_BMA.a_alpha = m*p_BMA.b_sigma; + p_BMA.a_alpha = m*p_BMA.b_alpha; else disp('VBA_MBA: Warning: mixture of deterministic and stochastic models!') p_BMA.b_alpha = Inf; diff --git a/stats&plots/mediationAnalysis0.m b/stats&plots/mediationAnalysis0.m index aa0e3009..16a0406b 100644 --- a/stats&plots/mediationAnalysis0.m +++ b/stats&plots/mediationAnalysis0.m @@ -50,7 +50,7 @@ [p1,stat1,df1,a1] = GLM_contrast([ones(n,1),X],Y,[0;1],'F',0); % regress M on X -[p2,stat2,df2,a2] = GLM_contrast([ones(n,1),M],Y,[0;1],'F',0); +[p2,stat2,df2,a2] = GLM_contrast([ones(n,1),X],M,[0;1],'F',0); % regress Y on both X and M [p3,stat3,df3,a3] = GLM_contrast([ones(n,1),X,M],Y,[0;1;0],'F',0); diff --git a/subfunctions/demo_CaBBI_FHN.m b/subfunctions/demo_CaBBI_FHN.m new file mode 100644 index 00000000..e795fd97 --- /dev/null +++ b/subfunctions/demo_CaBBI_FHN.m @@ -0,0 +1,253 @@ +% This script inverts the FHN model for the given fluorescence trace using +% CaBBI method: [ Rahmati. et al. 2016, PLoS Comput Biol 12(2) ] +% +% To user: +% I) in Step 1, determine the directory to the fluorescence trace (a real row vector); +% as default the code will download and use our published data for CaBBI +% II) in Step 2, determine the sampling rate (Hz) used for recording the fluorescence trace +% III) run the code +% IV) the main outputs of the CaBBI are listed in Step 11. +% +% Further options: +% 1) In Step 10, if necessary, you can change the spike-detection threshold which is used to extract +% spike-times from inferred voltage traces. Such change may increase the spike-detection accuracy. +% 2) In Step 6, if necessary, you can change some internal parameters of the optimization process. +% 3) For further description of the toolbox options, please see: +% http://mbb-team.github.io/VBA-toolbox/wiki/Controlling-the-inversion-using-VBA-options/ +% http://mbb-team.github.io/VBA-toolbox/wiki/VBA-output-structure/ + + +clear all; close all + +%% Step 1. Importing the fluorescence trace +% import the fluorescence trace +dataFolder = [pwd filesep 'Data' filesep]; +addpath('dataFolder'); +Fluor_trace_name = 'fluorescence_data8'; % file name of the fluorescence trace +Fluor_trace_dir = [dataFolder,Fluor_trace_name,'.mat']; + +if exist(Fluor_trace_dir,'file') == 0 % check whether the fluorescence data have been already downloaded + url = 'https://github.com/VahidRahmati/Sample_Data/blob/master/Sample_Data.zip?raw=true'; + urlwrite(url,[pwd filesep 'Data.zip']); % download the fluorescence data + unzip('Data.zip',[pwd filesep 'Data']); +end + +ydata = importdata(Fluor_trace_dir); % the fluorescence trace before pre-processing (a real row vector) + + +%% Step 2. Determining the sampling rate used to record the fluorescence trace +sampling_rate = 22.6; % in [Hz], for our data enter 22.6 +display(['the sampling rate used to record the fluorescence trace is << ',num2str(sampling_rate),' Hz >>']); + + +%% Step 3. Pre-processing +% prparing the data +warning('on') +if ~all(isfinite(ydata)) % check whether there is any NaN and/or Inf in the given fluorescence trace + ydata(~isfinite(ydata)) = []; % reovme the NaN and/or Inf elements from the trace + warning('the given fluorescence traces had some NaN and/or Inf elements, which were removed') +end + +if ~isreal(ydata) % check whether the given fluorescence trace contains any complex number + error('..... the given fluorescence traces contains complex numbers .....') +end + +% normalizing +if size(ydata,1) ~= 1; ydata = ydata'; end % fluorescence trace should be a row vector +ydata = (ydata - mean(ydata))./(max(mean(ydata),1)); % relativer fluorescncee trace +ydata = ydata - median(ydata); % removing the median + +% removing the slow drifts +xdata = 1:length(ydata); % vector of frames +degree = '4'; % Polynomial degree (4 or 3) +ft = fittype( ['poly',degree] ); % define the model to fit the slow drifts +opts = fitoptions( ft ); +opts.Lower = [-Inf -Inf -Inf -Inf -Inf]; +opts.Robust = 'Bisquare'; +opts.Upper = [Inf Inf Inf Inf Inf]; +warning('off') +[fitresult, gof] = fit( xdata', ydata', ft, opts ); % fitting the polynomial model to fluorecence trace +warning('on') +Coeffs = coeffvalues(fitresult); +if degree == '4' + fittedcurve = Coeffs(1)*xdata.^4 + Coeffs(2)*xdata.^3 + Coeffs(3)*xdata.^2 + Coeffs(4)*xdata + Coeffs(5); +elseif degree == '3' + fittedcurve = Coeffs(1)*xdata.^3 + Coeffs(2)*xdata.^2 + Coeffs(3)*xdata.^1 + Coeffs(4); +end + +Fluor_trace = ydata - fittedcurve; % pre-processed fluorescence trace + + +%% Step 4. Scaling the data +scale = 1; % preferably set it a number between e.g. 0.4 and 1 +Fluor_trace = scale* Fluor_trace/max(Fluor_trace); +plot(xdata, scale*ydata/max(ydata) ,'b',xdata, Fluor_trace(1:numel(Fluor_trace)),'r'); +legend('original','pre-processed') +xlabel('time [frame]') +ylabel('a.u.') +title('fluorescence trace') + + +%% Step 6. The generative model +% assigning the generative model +f_fname = @f_CaBBI_FHN; % evolution function +g_fname = @g_CaBBI; % observation function + +% setting the parameters of the FHN model used for inversion +nFrames = numel(Fluor_trace); % number of recording frames +dt = 0.2; % [ms], time step size of the FHN model +inF.dt = dt; +inF.k_Ca0 = 0.002; % scale parameter of [Ca2+] kinetics +TauCa_realTime = 7500; % [ms], the prior for decay time-constant of calcium transients in real time +dt_real = 1000/sampling_rate; % [ms], frame duration (temporal precision) of fluorescence data +inF.Tau_Ca0 = TauCa_realTime*dt/dt_real; % effective prior mean of tau_Ca in inversion; see Eqn. 17 +inG.k_F0 = 5; % scale parameter of the fluorescence observatios +inG.ind = 3; % index of [Ca2+] variable + + +%% Step 6. Setting the internal parameters of the optimization (VB-Laplace) method +options.backwardLag = 2; % lag parameter (selcet either 2 or 3); increasing it e.g. to 3 (or higher) "may" provide more accurate inversions +options.updateHP = 1; % this allows the precison (hyper) parameters to be learnt from data +options.MaxIterInit = 0; +options.MaxIter = 2; % or e.g. 3 % maximum number of optimization iterations; increasing it will lead to a + % complete optimization process, although it may take long time to be converged) + +% options.GnFigs = 0; % set to zero to speed up the inversion +options.DisplayWin = 0; % set to zero to speed up the inversion +options.verbose = 1; % set to zero to speed up the inversion +options.decim = 1; % this determine the micro-time resolution; increasing it e.g. to 2 (or higher) "may" provide more accurate inversions + + +%% Step 7. Assigning the prior densitities: (mu --> mean), (Sigma --> variance) +% prior densities for initial conditons +priors.muX0 = [-1.2; -0.62; 0.05]; % [V(0), W(0), [Ca2+](0)] +priors.SigmaX0(1,1) = 4; +priors.SigmaX0(2,2) = 1; +priors.SigmaX0(3,3) = 25; + +% prior densities for evolution parameters +priors.muTheta = [0; 0; 0]; +priors.SigmaTheta(1,1) = 5; % for X_Tau_Ca +priors.SigmaTheta(2,2) = 100; % for Ca_base +priors.SigmaTheta(3,3) = 5; % for X_k_Ca + +% prior densities for observation parameters +priors.muPhi = [min(Fluor_trace)/2;0]; +priors.SigmaPhi(1,1) = 0.25; % for d_F +priors.SigmaPhi(2,2) = 1; % for X_k_F + +% prior densities for precision parameters: Gamma distributions +priors.a_alpha = 1e0; +priors.b_alpha = 1e0; +priors.a_sigma = 1e0; +priors.b_sigma = 1e0; + +% building the precison matrix for evolution variables +priors.iQx = cell(nFrames,1); +for t = 1:nFrames + dq(1) = (2); % on V + dq(2) = (30); % on W + dq(3) = 10; % on [Ca2+] + priors.iQx{t} = diag(dq); +end + +% building the precision verctor the fluorescence observations +priors.iQy = cell(nFrames,1); +for t=1:nFrames + priors.iQy{t} = 3; % on F +end + +% setting the option structure +options.inF = inF; +options.inG = inG; +options.priors = priors; + +% assigning the dimensions for parameter space +dim.n = 3; % number of neuroal variables +dim.n_phi = 2; % number of observation paramters +dim.n_theta = 3; % number of evolution parameters + + +%% Step 8. Inversion: inverting the FHN model for the given fluorescence trace +[posterior, out] = VBA_NLStateSpaceModel(Fluor_trace,zeros(1,numel(Fluor_trace)),f_fname,g_fname,dim,options); + + +%% Post-processing +%% Step 9. Extracting the estimated Voltage and [Ca2+] traces (i.e. the posterior means) + +V_inferred = posterior.muX(1,1:end); % inferred Voltage (posterior mean) +W_inferred = posterior.muX(2,1:end); % inferred recovery variable (posterior mean) +Ca_inferred = posterior.muX(end,1:end); % inferred [Ca2+] kinetics (posterior mean) + +% Plotting the estimated voltage and [Ca2+] traces +Calcium_basal = 50; % in [nM], the basal [Ca2+] concentration + +% Plot the pre-processed fluorescence trace +figure +subplot(411); plot(1:nFrames, Fluor_trace,'g','LineWidth',2) +legend('fluorescence') +set(gca,'xLim',[0 nFrames]) +axis tight + +% Plot the posterior mean of [Ca2+] trace +subplot(412); plot(1:nFrames, Ca_inferred+Calcium_basal,'r','LineWidth',2) +legend('inferred [Ca2+]') +set(gca,'xLim',[0 nFrames]) +axis tight + +% plot the posterior mean of voltage trace +subplot(413); plot(1:nFrames, V_inferred,'b','LineWidth',2) +set(gca,'xLim',[0 nFrames]) +axis tight + + +%% Step 10. Detecting spike times (or event onsets) from inferred voltage trace +Detection_threshold = 0; % the spike (or event) detection threshold + % (if necessary, change the detection threshold to e.g. -0.2 or 0.5) +idxSpikes = zeros(1,numel(V_inferred)); +count = 1; +for j=1:numel(V_inferred)-1 + + if( V_inferred(j+1)>Detection_threshold && V_inferred(j) 1) && ( (j-indices(count-1)) > ceil(6/dt) ) + idxSpikes(j) = j; + indices(count) = j; + count = count+1; + end + + if count == 1 + idxSpikes(j) = j; + indices(count) = j; + count = count+1; + end + + end + +end +idxSpikeToShow = idxSpikes; +idxSpikeToShow(idxSpikeToShow~=0) = 1; + +subplot(414); plot(1:nFrames,idxSpikeToShow,'b','LineWidth',2) +legend('inferred spiketimes'); +xlabel('time [frame]') +set(gca,'xLim',[0 nFrames]) +axis tight + +subplot(413); hold on; plot(1:nFrames, Detection_threshold*ones(1,nFrames),'-.k'); +legend('inferred voltage','threshold'); hold off + +% extracted spike (or event) times +SpikeTimes_Frame = idxSpikeToShow(idxSpikeToShow~=0)'; % reconstructed spike (or event) times, in [frame] +SpikeTimes_millisecond = SpikeTimes_Frame*1000/sampling_rate; % reconstructed spike (or event) times, in [ms] + + +%% Step 11. Main outputs of CaBBI method +CaBBI.nFrames = nFrames; % the number of recordings frames of the given fluorescence trace +CaBBI.sampling_rate = sampling_rate; % the sampling rate used for recording the fluorescence trace, in [Hz] +CaBBI.Fluor_trace = Fluor_trace; % the pre-processed fluorescebce trace +CaBBI.V_inferred = V_inferred; % inferred Voltage trace (posterior mean) +CaBBI.W_inferred = W_inferred; % inferred trace of recovery variable (posterior mean) +CaBBI.Ca_inferred = Ca_inferred; % inferred [Ca2+] trace (posterior mean) +CaBBI.SpikeTimes_Frame = SpikeTimes_Frame; % reconstructed spike (or event) times, in [frame] +CaBBI.SpikeTimes_millisecond = SpikeTimes_millisecond; % reconstructed spike (or event) times, in [ms] diff --git a/subfunctions/demo_CaBBI_QGIF.m b/subfunctions/demo_CaBBI_QGIF.m new file mode 100644 index 00000000..8c036ca1 --- /dev/null +++ b/subfunctions/demo_CaBBI_QGIF.m @@ -0,0 +1,250 @@ +% This script inverts the QGIF model for the given fluorescence trace using +% CaBBI method: [ Rahmati. et al. 2016, PLoS Comput Biol 12(2) ] +% +% To user: +% I) in Step 1, determine the directory to the fluorescence trace (a real row vector); +% as default the code will download and use our published data for CaBBI +% II) in Step 2, determine the sampling rate (Hz) used for recording the fluorescence trace +% III) run the code +% IV) the main outputs of the CaBBI are listed in Step 11. +% +% Further options: +% 1) In Step 10, if necessary, you can change the spike-detection threshold which is used to extract +% spike-times from inferred voltage traces. Such change may increase the spike-detection accuracy. +% 2) In Step 6, if necessary, you can change some internal parameters of the optimization process. +% 3) For further description of the toolbox options, please see: +% http://mbb-team.github.io/VBA-toolbox/wiki/Controlling-the-inversion-using-VBA-options/ +% http://mbb-team.github.io/VBA-toolbox/wiki/VBA-output-structure/ + + +clear all; close all + +%% Step 1. Importing the fluorescence trace +% import the fluorescence trace +dataFolder = [pwd filesep 'Data' filesep]; +addpath('dataFolder'); +Fluor_trace_name = 'fluorescence_data8'; % file name of the fluorescence trace +Fluor_trace_dir = [dataFolder,Fluor_trace_name,'.mat']; + +if exist(Fluor_trace_dir,'file') == 0 % check whether the fluorescence data have been already downloaded + url = 'https://github.com/VahidRahmati/Sample_Data/blob/master/Sample_Data.zip?raw=true'; + urlwrite(url,[pwd filesep 'Data.zip']); % download the fluorescence data + unzip('Data.zip',[pwd filesep 'Data']); +end + +ydata = importdata(Fluor_trace_dir); % the fluorescence trace before pre-processing (a real row vector) + + +%% Step 2. Determining the sampling rate used to record the fluorescence trace +sampling_rate = 22.6; % in [Hz], for our data enter 22.6 +display(['the sampling rate used to record the fluorescence trace is << ',num2str(sampling_rate),' Hz >>']); + + +%% Step 3. Pre-processing +% prparing the data +warning('on') +if ~all(isfinite(ydata)) % check whether there is any NaN and/or Inf in the given fluorescence trace + ydata(~isfinite(ydata)) = []; % reovme the NaN and/or Inf elements from the trace + warning('the given fluorescence traces had some NaN and/or Inf elements, which were removed') +end + +if ~isreal(ydata) % check whether the given fluorescence trace contains any complex number + error('..... the given fluorescence traces contains complex numbers .....') +end + +% normalizing +if size(ydata,1) ~= 1; ydata = ydata'; end % fluorescence trace should be a row vector +ydata = (ydata - mean(ydata))./(max(mean(ydata),1)); % relativer fluorescncee trace +ydata = ydata - median(ydata); % removing the median + +% removing the slow drifts +xdata = 1:length(ydata); % vector of frames +degree = '4'; % Polynomial degree (4 or 3) +ft = fittype( ['poly',degree] ); % define the model to fit the slow drifts +opts = fitoptions( ft ); +opts.Lower = [-Inf -Inf -Inf -Inf -Inf]; +opts.Robust = 'Bisquare'; +opts.Upper = [Inf Inf Inf Inf Inf]; +warning('off') +[fitresult, gof] = fit( xdata', ydata', ft, opts ); % fitting the polynomial model to fluorecence trace +warning('on') +Coeffs = coeffvalues(fitresult); +if degree == '4' + fittedcurve = Coeffs(1)*xdata.^4 + Coeffs(2)*xdata.^3 + Coeffs(3)*xdata.^2 + Coeffs(4)*xdata + Coeffs(5); +elseif degree == '3' + fittedcurve = Coeffs(1)*xdata.^3 + Coeffs(2)*xdata.^2 + Coeffs(3)*xdata.^1 + Coeffs(4); +end + +Fluor_trace = ydata - fittedcurve; % pre-processed fluorescence trace + + +%% Step 4. Scaling the data +scale = 1; % preferably set it a number between e.g. 0.4 and 1 +Fluor_trace = scale* Fluor_trace/max(Fluor_trace); +plot(xdata, scale*ydata/max(ydata) ,'b',xdata, Fluor_trace(1:numel(Fluor_trace)),'r'); +legend('original','pre-processed') +xlabel('time [frame]') +ylabel('a.u.') +title('fluorescence trace') + + +%% Step 5. The generative model +% assigning the generative model +f_fname = @f_CaBBI_QGIF; % evolution function +g_fname = @g_CaBBI; % observation function + +% setting the parameters of the QGIF model used for inversion +nFrames = numel(Fluor_trace); % number of recording frames +dt = 0.2; % [ms], time-step size of the QGIF model +inF.dt = dt; +inF.k_Ca0 = 0.002; % scale parameter of [Ca2+] kinetics +TauCa_realTime = 7500; % [ms], the prior for decay time-constant of calcium transients in real time +dt_real = 1000/sampling_rate; % [ms], frame duration (temporal precision) of fluorescence data +inF.Tau_Ca0 = TauCa_realTime*dt/dt_real; % efective prior mean for tau_Ca in inversion; see Eqn. 17 +inG.k_F0 = 5; % scale parameter of the fluorescence observations +inG.ind = 2; % index of [Ca2+] variable + + +%% Step 6. Setting the internal parameters of the optimization (VB-Laplace) method +options.backwardLag = 2; % lag parameter (selcet either 2 or 3); increasing it e.g. to 3 (or higher) "may" provide more accurate inversions +options.updateHP = 1; % this allows the precison (hyper) parameters to be learnt from data +options.MaxIterInit = 0; +options.MaxIter = 2; % or e.g. 3 % maximum number of optimization iterations; increasing it will lead to a + % complete optimization process, although it may take long time to be converged) + +% options.GnFigs = 0; % set to zero to speed up the inversion +options.DisplayWin = 0; % set to zero to speed up the inversion +options.verbose = 1; % set to zero to speed up the inversion +options.decim = 1; % this determine the micro-time resolution; increasing it e.g. to 2 (or higher) "may" provide more accurate inversions + + +%% Step 7. Assigning the prior densities: (mu --> mean), (Sigma --> variance) +% prior densities for initial conditions +priors.muX0 = [ -63.62; 5.65]; % [V(0), [Ca2+](0)] +priors.SigmaX0(1,1) = 100; +priors.SigmaX0(2,2) = 25; + +% prior densities for evolution parameters +priors.muTheta = [0; 0; 0]; +priors.SigmaTheta(1,1) = 5; % for X_Tau_Ca +priors.SigmaTheta(2,2) = 100; % for Ca_base +priors.SigmaTheta(3,3) = 5; % for X_k_Ca + +% prior densities for observation parameters +priors.muPhi = [min(Fluor_trace)/2; 0]; +priors.SigmaPhi(1,1) = 0.25; % for d_F +priors.SigmaPhi(2,2) = 1; % for X_k_F + +% prior densities for precision parameters: Gamma distributions +priors.a_alpha = 1e0; +priors.b_alpha = 1e0; +priors.a_sigma = 1e0; +priors.b_sigma = 1e0; + +% building the precision matrix for evolution variables + priors.iQx = cell(nFrames,1); +for t = 1:nFrames + dq(1) = 5e-3; % on V + dq(2) = 10; % on [Ca2+] + priors.iQx{t} = diag(dq); +end + +% building the precision vector for fluorescence observations +priors.iQy = cell(nFrames,1); +for t=1:nFrames + priors.iQy{t} = (3); % on F +end + +% setting the option structure +options.inF = inF; +options.inG = inG; +options.priors = priors; + +% assigning the dimensions for parameter space +dim.n = 2; % number of neural variables +dim.n_phi = 2; % number of observation parameters +dim.n_theta = 3; % number of evolution parameters + + +%% Step 8. Inversion: inverting the QGIF model for the given fluorescence trace +[posterior, out] = VBA_NLStateSpaceModel(Fluor_trace,zeros(1,numel(Fluor_trace)),f_fname,g_fname,dim,options); + + +%% Post-processing +%% Step 9. Extracting the estimated Voltage and [Ca2+] traces (i.e. the posterior means) + +V_inferred = posterior.muX(1,1:end); % inferred Voltage (posterior mean) +Ca_inferred = posterior.muX(end,1:end); % inferred [Ca2+] kinetics (posterior mean) + +% Plotting the estimated voltage and [Ca2+] traces +Calcium_basal = 50; % in [nM], the basal [Ca2+] concentration + +% Plot the pre-processed fluorescence trace +figure +subplot(411); plot(1:nFrames, Fluor_trace,'g','LineWidth',2) +legend('fluorescence') +set(gca,'xLim',[0 nFrames]) +axis tight + +% Plot the posterior mean of [Ca2+] trace +subplot(412); plot(1:nFrames, Ca_inferred+Calcium_basal,'r','LineWidth',2) +legend('inferred [Ca2+]') +set(gca,'xLim',[0 nFrames]) +axis tight + +% plot the posterior mean of voltage trace +subplot(413); plot(1:nFrames, V_inferred,'b','LineWidth',2) +set(gca,'xLim',[0 nFrames]) +axis tight + + +%% Step 10. Detecting spike times (or event onsets) from inferred voltage trace +Detection_threshold = 0; % the spike (or event) detection threshold + % (if necessary, change the detection threshold to e.g. -10) +idxSpikes = zeros(1,numel(V_inferred)); +count = 1; +for j=1:numel(V_inferred)-1 + + if( V_inferred(j+1)>Detection_threshold && V_inferred(j) 1) && ( (j-indices(count-1)) > ceil(6/dt) ) + idxSpikes(j) = j; + indices(count) = j; + count = count+1; + end + + if count == 1 + idxSpikes(j) = j; + indices(count) = j; + count = count+1; + end + + end + +end +idxSpikeToShow = idxSpikes; +idxSpikeToShow(idxSpikeToShow~=0) = 1; + +subplot(414); plot(1:nFrames,idxSpikeToShow,'b','LineWidth',2) +legend('inferred spiketimes'); +xlabel('time [frame]') +set(gca,'xLim',[0 nFrames]) +axis tight + +subplot(413); hold on; plot(1:nFrames, Detection_threshold*ones(1,nFrames),'-.k'); +legend('inferred voltage','threshold'); hold off + +% extracted spike (or event) times +SpikeTimes_Frame = idxSpikeToShow(idxSpikeToShow~=0)'; % reconstructed spike (or event) times, in [frame] +SpikeTimes_millisecond = SpikeTimes_Frame*1000/sampling_rate; % reconstructed spike (or event) times, in [ms] + + +%% Step 11. Main outputs of CaBBI method +CaBBI.nFrames = nFrames; % the number of recordings frames of the given fluorescence trace +CaBBI.sampling_rate = sampling_rate; % the sampling rate used for recording the fluorescence trace, in [Hz] +CaBBI.Fluor_trace = Fluor_trace; % the pre-processed fluorescebce trace +CaBBI.V_inferred = V_inferred; % inferred Voltage trace (posterior mean) +CaBBI.Ca_inferred = Ca_inferred; % inferred [Ca2+] trace (posterior mean) +CaBBI.SpikeTimes_Frame = SpikeTimes_Frame; % reconstructed spike (or event) times, in [frame] +CaBBI.SpikeTimes_millisecond = SpikeTimes_millisecond; % reconstructed spike (or event) times, in [ms] + diff --git a/subfunctions/f_CaBBI_FHN.m b/subfunctions/f_CaBBI_FHN.m new file mode 100644 index 00000000..24685772 --- /dev/null +++ b/subfunctions/f_CaBBI_FHN.m @@ -0,0 +1,44 @@ +% FHN evolution function: +% +% [fx] = f_CaBBI_FHN(Xt,Theta,Isyn,inF) +% +% this function contains the equations and parameters of the FHN neuron +% model, and the over-time integration scheme + + +function [fx] = f_CaBBI_FHN(Xt,Theta,Isyn,inF) + + +dt = inF.dt; +Tau_Ca0 = inF.Tau_Ca0; + + +% Hidden variables +V = Xt(1); % voltage +W = Xt(2); % recovery variable +Ca = Xt(3); % [Ca2+] kinetics + +% evolution parameters +Tau_Ca = Tau_Ca0* exp(Theta(1)); % decay time-constant of [Ca2+] kinetics +Ca_base = Theta(2); % basal value of [Ca2+] +scale = inF.k_Ca0* exp(Theta(3)); % scale parameter of [Ca2+] kinetis + +% parameters of the evolution mode +g_Ca = 5; % maximal conductance of the calcium channel +E_Ca = 120; % reversal potential of the calcium channel + +% compute the current of HVA calcium channel +V_scaled = (30*V) - 40; +s_inf = 1./( 1 + exp(- (V_scaled+45)/5 ) ); +I_Ca = g_Ca .* s_inf.* (V_scaled-E_Ca); % HVA calcium current + + +% build the right hand side (RHS) of the evolution's differential equations +xdot = [ ( V - (power(V,3)/3) - W + Isyn ) + (0.08*(V+0.7-(0.8*W))) + ( -scale*I_Ca - ( (Ca-Ca_base)/Tau_Ca) ) + ]; + + +fx = Xt + dt.*xdot; + diff --git a/subfunctions/f_CaBBI_QGIF.m b/subfunctions/f_CaBBI_QGIF.m new file mode 100644 index 00000000..096767cf --- /dev/null +++ b/subfunctions/f_CaBBI_QGIF.m @@ -0,0 +1,56 @@ +% QGIF evolution function: +% +% [fx] = f_CaBBI_QGIF(Xt,Theta,Isyn,inF) +% +% this function contains the equations and parameters of the QGIF neuron +% model, and the over-time integration scheme + + +function [fx] = f_CaBBI_QGIF(Xt,Theta,Isyn,inF) + +dt = inF.dt; +Tau_Ca0 = inF.Tau_Ca0; + +% hidden variables +V = Xt(1); +Ca = Xt(2); + +% evolution parameters +Ca_base = Theta(2); % basal value of [Ca2+] +scale = inF.k_Ca0*exp(Theta(3)); % scale parameter of [Ca2+] kinetics +Tau_Ca = Tau_Ca0* exp(Theta(1)); % decay time-constant of [Ca2+] kinetics + + +% other parameters of the evolution model +Cm = 1; % membrane capacitance +g_L = 0.05; % leak conductance +g_Ca = 5; % maximal conductance of HVA calcium channel +E_Ca = 120; % reversal potential of HVA calcium channel +E_L = -65; % reversal potential of leak current +Vth = -59.9; % voltage threshold +Ith = 0.16; % threshold current +deltaV = 3.48; % spike slope factor + + +% Compute the currents +s_inf = 1./( 1 + exp(- (V+45)/5 ) ); +I_Ca = g_Ca .* s_inf.* (V-E_Ca); % HVA calcium current +I_L = g_L * (V-E_L); % leak current +f_sp = ((g_L/(2*deltaV))*power(V-Vth,2) ) + I_L -Ith; +V_peak = 30; +I_rep = 120; +b_rep = 1e8; +c_rep = 0.1; +a_rep = 1e12; +sigma_peak = 1; +G_rep = a_rep*(1/(sqrt(pi)*sigma_peak))*exp(- power(V - V_peak,2)/power(sigma_peak,2) ); +f_rep =I_rep* 1./ (1+exp(-b_rep*((G_rep)-c_rep))); + + +% build the right hand side (RHS) of the evolution's differential equations +xdot = [ ((-I_L+ f_sp + Isyn) ./ Cm - (f_rep/dt) ) + ( -scale*I_Ca - ( (Ca-Ca_base)/Tau_Ca) ) + ]; + + +fx = Xt + dt.*xdot; diff --git a/subfunctions/g_CaBBI.m b/subfunctions/g_CaBBI.m new file mode 100644 index 00000000..bd071ba3 --- /dev/null +++ b/subfunctions/g_CaBBI.m @@ -0,0 +1,18 @@ +% The observation function: +% +% function [F] = g_CaBBI(Xt,Phi,I_inp,inG) +% +% this function maps [Ca2+] kinetics to fluorescence observations (F) +% through a non-linear saturating function + + +function [F] = g_CaBBI(Xt,Phi,I_inp,inG) + +ind = inG.ind; +Ca = Xt(ind); % [Ca2+] kinetics +d_F = Phi(1); % offset parameter +Kd = 200; % dissociation constant +scale = inG.k_F0*exp(Phi(2)); % scale paraemeter +F = scale*(Ca/(Ca + Kd)) + d_F; % nonlinear mapping from [Ca2+] to F + +end diff --git a/subfunctions/getHyperpriors.m b/subfunctions/getHyperpriors.m deleted file mode 100644 index 09f94531..00000000 --- a/subfunctions/getHyperpriors.m +++ /dev/null @@ -1,40 +0,0 @@ -function [a, b]=getHyperpriors(variance,p_min,p_max) -%% compute hyperpriors parameters given you want to eplain between p_min -% and p_max fraction of the total variance. -%% ex. -% [a, b]=getHyperpriors(200,.1,.9); -% figure; -% r=random('gam',a,1/b,1,5000); -% hist(r,50); -% hold on -% plot([prec_min, prec_max],[100 100],'r'); - -%% check parameters -if nargin==1 - p_min=.1; - p_max=.9; -elseif nargin==3 - if any([p_min, p_max]<0) || any([p_min, p_max]>1) - error('*** getHyperpriors: p_min and p_max should be between 0 and 1.'); - end - if p_min > p_max - error('*** getHyperpriors: p_min should be inferior to p_max.'); - end -else - error('*** getHyperpriors: wrong number of argument'); -end - -p_max = min(p_max, 1-eps); - -% residual variance -var_min = variance*(1-p_max); -var_max = variance*(1-p_min); - -% expressed as precision -prec_min = 1/var_max; -prec_max = 1/var_min; - -% approximate 98% confidence interval -a = 6*((prec_min+prec_max)^2)/((prec_max-prec_min)^2); -b = 12*(prec_min+prec_max)/((prec_max-prec_min)^2); -