NMDA and CMM_NMDA models are not initialised at steady state in DCM for ERP & CSD #121

The neural states for two models: NMDA and CMM_NMDA are not initialised at their steady state before the model sees any input. As a result, these models start from arbitrary initial conditions which affects downstream inversion.

Root cause and call chain:

1. spm_dcm_erp:L222 and spm_dcm_csd:L117 call spm_dcm_x_neural to initialise the specified model.
2. For CMM_NMDA, spm_dcm_x_neural calls spm_x_cmm_NMDA(P), and for NMDA, it calls spm_x_nmda(P)
3. Both spm_x_cmm_NMDA and spm_x_nmda do some prep and call spm_dcm_neural_x(P,M) which is the steady-state solver
4. spm_dcm_neural_x has cases defined for CMM and NMM (via MFM) but none for CMM_NMDA and NMDA
5. As a result, no Newton-Raphson iterations take place for those two models and the function returns the naive initial conditions, which for CMM_NMDA correspond to -50mV membrane potential and 1/8 conductance, and for NMDA, -70mV and 1/8 respectively.

Illustrative minimal reproducible example

The following code snippet loops over 6 models and simulates evoked LFP for a single node under default conditions

% Loop over these models
models = {'ERP', 'NMM', 'NMDA', 'CMC', 'CMM', 'CMM_NMDA'};
nmodels = numel(models);

f = figure('Position',[0,0,1600,600]);
for m=1:nmodels

    % Initialize a single node
    A = {0, 0, 0}; B = {0}; C = 1;
    Ns = numel(diag(A{1}));
    model = models{m};
    
    % Get default priors, and intialise model using them
    pE = spm_dcm_neural_priors(A,B,C,model);
    [x,f,h] = spm_dcm_x_neural(pE,model);
    
    % Use priors as parameter set for simulation
    P = pE;
    
    % Initialise model struct with minimum necessary variables
    M = [];
    M.model     = model;                
    M.x         = x;                    % Initial states
    M.f         = f;                    % 'fx' function file
    M.m         = numel(C);             % Number of nodes
    M.n         = length(spm_vec(M.x)); % Number of states
    M.ons       = 60;                   % Onset of input (ms)
    M.dur       = 16;                   % Duration of input (ms, Gaussian SD)
    M.ns        = 500;                  % Number of samples
    
    % Input
    U = []; U.dt = 0.001;               % Sampling interval (s)
    pst = 1:M.ns;                       % Peri-stimulus window
    
    % Generate ERP
    yy = spm_gen_erp(P,M,U);
    
    % Initialize a dummy dipfit struct
    dipfit = []; dipfit.model    = model;
    dipfit.type = 'LFP'; dipfit.location = 0; dipfit.Ns = 1; dipfit.Nc = 1;
    
    % Prepare 'LFP' output of model
    gE = spm_L_priors(dipfit); % Get observation model parameters
    L = spm_lx_erp(gE, dipfit); % Obtain leadfield (same as gE.J in this case)
    x0  = ones(Ns,1)*spm_vec(M.x)'; % Get expansion points for all states
    K = yy{1} - x0; % Center generated states on expansion point
    y = K*L'; % Project on to lead field to get LFP output of entire node
    
    % Plot
    subplot(2,nmodels,m)
    hold on; plot(pst,y); yline(0); hold off
    xlim([pst(1), pst(end)])
    title(['Evoked LFP for ' spm_str_manip(model, 'x')])
    xlabel('Time (ms)')
    if m==1; ylabel('Activity (AU)'); end
    subplot(2,nmodels,nmodels+m)
    bar(x0)
    title(['Initial states for ' spm_str_manip(model, 'x')])
end

The plot produced by this (on spm main as of 5fa3fb2) looks like:

As is evident in the top row, the ERP, NMM, CMC and CMM models show a biphasic response that returns to 0, whereas NMDA does show a biphasic response which returns to a non-zero baseline and CMM_NMDA does neither. The bottom row shows the initial states for each model (M.x) which are specified for the conductance-based models (hence 0 for all states of ERP and CMC). The initial states for NMM and CMM are slightly tweaked from the baseline values specified for them, while those for NMDA and CMM_NMDA have not.

Fixes (Part 1)

The most obvious fix involves just adding the necessary case arguments in spm_dcm_neural_x for those two models so that Netwon-Raphson can work. That works for CMM_NMDA but not for NMDA (and will crash) because certain fields are not in the right shape (M.x is a cell), and critically, the Jacobian computed in spm_fx_nmda:L276 is incorrectly obtained by spm_diff on spm_fx_mfm instead of spm_fx_nmda

The three fixes then are:

1. add spm_fx_cmm_NMDA and spm_fx_nmda to spm_dcm_neural_x
2. handle M.x in spm_x_nmda, and
3. compute Jacobian correctly in spm_fx_nmda.

Applying these and repeating the snippet above now produces this plot:

Now NMDA and CMM_NMDA are both returning to 0 (top row) and their initial states are different from their default values. But CMM_NMDA is still not producing a biphasic response.

Fixes (Part 2)

Comparing spm_fx_cmm_NMDA and spm_fx_cmm reveals many differences beyond the inclusion of the NMDA receptor, especially in intrinsic connectivity between the 4 populations, their membrane capacitance and rate constants for receptors. They are fundamentally different models, but one difference that stands out is that the input is scaled by 1000 for CMM but not for CMM_NMDA. This might be the pathology, and scaling the input similarly for CMM_NMDA appears to fix it. Here's the plot from snippet above:

Now CMM_NMDA also produces a biphasic response that returns to 0.

Open issues

• With the above fixes, the peak amplitudes of ERP, CMC, CMM and CMM_NMDA are in a similar range, but those of NMM and NMDA are quite large. In their respective spm_fx* functions, inputs are not scaled at all for NMDA while they are scaled by 2 for NMM, which may not be enough. Should their inputs be scaled as well? If so, to what extent. (For e.g. scaling NMDA's input by 1000 does produce a more realistic-looking waveform with weaker damping, but results in very low amplitude output for NMM )