Fix issues affecting Deco2018 Figure 3 reproduction

[zhouzx-UCAS]

This PR fixes two implementation issues identified during a detailed reproduction of Deco et al., 2018 (Current Biology), especially Figure 3B.

The main fix clips the coupling state to [0, 1] before applying linear coupling. This restores consistency with the WholeBrain implementation, where Se/Si are clamped before long-range coupling is computed. Without this step, neuronumba produces systematically different dynamics and cannot reproduce Figure 3B correctly.

In addition, this PR fixes several issues in the Deco2014 FIC implementation uncovered during the reproduction process.

[neich]

Hi,

After looking at your PR I have two comments.

First, you are doing the clipping at LinearModel, which I am not sure it is correct, since it is the general coupling method for all linear models, and you do not now at that point what is the nature of the state variables. You do not know if they need to be clipped. In a different model than Deco2014 maybe they have another range. Also, model Deco2014 already clips the gating variables. Look at this

And second, wrt to fic mechanism implemented in FICDeco2014, I am sure your changes are right. But I strongly not recommend to use FICDeco2014 as fic mechanism since it is a poorly implemented optimization algorithm, and more important you have FICHerzog2022 mechanism that computes exactly the same, even more accurately, and extremely faster. I have to add a deprecated information on FICDeco2014. In fact, all models that use automatic fic use FICHerzog2022.

Let me know what you think @zhouzx-UCAS

[neich]

@zhouzx-UCAS Also look at the latest PRs from 3 weeks ago, I fixed some of the neuronumba previuos problems, and the Deco 2018 figures are ok now. I have checked your fork and your main branch is 49 commits behind the neuronumba main branch.

[zhouzx-UCAS]

Hi Ignacio,

Thank you for the quick and thoughtful reply! I agree with both of your comments, but I still have a few thoughts that I would like to discuss with you.

0. Clarification about the code version used for testing

I would like to clarify that, over the past few weeks, I have been running all of my tests on the version that already includes your latest merged PRs. I know that the main branch of my fork is not up to date, but all my tests and pushes were done on the PR branch:

https://github.com/zhouzx-UCAS/neuronumba/tree/PR

I am aware that you have already made many important updates to the Deco 2018-related implementation, and I would not have been able to reach the current progress without them. For example, I noticed that the current latest version uses a noise level of 1e-2, and in my comparison tests, this also appears to be an important factor in allowing the FCD distance in Figure 3 to decrease substantially. Please let me know if I am mistaken about that.

1. About clipping before coupling

Regarding clipping before coupling, it does not seem necessary for reproducing Figure 3A, but it does seem necessary for Figure 3B, regardless of whether auto_fic is used. You mentioned that clipping at LinearModel may be inappropriate for some models, and I think that is a very important point that I had not considered before (as over the past period I had focused entirely on Deco2014 and its derived models). In that sense, I agree that the current PR is not the right form for this fix, and I will close it.

However, I did see that the current version clips the gating variables in the Deco2014/2018 models, and for that reason I did not suspect clipping to be the cause during nearly two weeks of investigation. The key point is that in the current implementation the gating variables are clipped inside dfun, but the coupling term is computed earlier and passed into dfun as an argument. So the current clipping does not constrain the state actually used for long-range coupling.

        def Deco2018_dfun(state: NDA_f8_2d, coupling: NDA_f8_2d):
            # Clamping, needed in Deco2018 model and derivatives...
            Se = state[0, :].clip(0.0,1.0)
            Si = state[1, :].clip(0.0,1.0)

            # Eq for I^E (5). I_external = 0 => resting state condition.
            Ie = m[np.intp(P.Jext_e)] * m[np.intp(P.I0)] + m[np.intp(P.w)] * m[np.intp(P.J_NMDA)] * Se + m[np.intp(P.J_NMDA)] * coupling[0, :] - m[np.intp(P.J)] * Si + m[np.intp(P.I_external)]

As you can see, the clipping there affects the local state variables used in the local terms, but it does not affect the value of the coupling term itself.

By contrast, in WholeBrain, the coupling is computed after clipping:

@jit(nopython=True)
def dfun(simVars, coupling, I_external):
    Se = simVars[0]; Si = simVars[1]  # should be [Se, Si] = simVars
    Se = utils.doClamping(Se); Si = utils.doClamping(Si)  # Clamping, needed in Deco2014 model and derivatives...

    coupl = coupling.couple(Se)
    Ie = I0 * Jexte + w * J_NMDA * Se + we * J_NMDA * coupl - J * Si + I_external  # Eq for I^E (5). I_external = 0 => resting state condition.

So although the difference looks small, the two implementations are not actually equivalent: in neuronumba, clipping currently happens too late to affect the coupling term, whereas in WholeBrain the long-range coupling is computed from the already clamped Se.

In practice, this difference seems to have only a limited effect in most cases, but it completely changes the Figure 3B pattern. I can also send you the comparison plots if you would like to check whether this difference is indeed fundamental.

If you agree, one possible solution would be to make coupling-state preprocessing model-specific: the base class would keep the default identity behavior, while models such as Deco2014/2018 could override a dedicated preprocessing hook to clip bounded gating variables before coupling is computed. If that sounds reasonable to you, I would be happy to submit a separate PR for that.

2. About FIC

Regarding FIC, I fully agree that FICHerzog2022 is the better mechanism and should be the default. I only fixed FICDeco2014 because, while trying to identify the source of the Figure 3B mismatch, I needed to control and compare each implementation step as strictly as possible. Without a working FICDeco2014, I could not determine where the pattern difference actually came from.

When using FICHerzog2022, the G value giving the best FCD in Figure 3A is slightly larger than 2.1 (around 2.4), and Figure 3B could not be reproduced before fixing the clipping-before-coupling issue. The good news is that I re-tested this today, and after fixing the clipping-before-coupling issue, Figure 3B can also be reproduced using FICHerzog2022 with the optimal G = 2.4. This suggests that the original Deco2018 finding is robust across different FIC mechanisms.

If you think it is still useful to keep FICDeco2014 available for users who want to accurately reproduce earlier results, I would also be happy to submit a separate PR for that, possibly with a deprecation note indicating that it is not recommended for general use.

@neich Let me know what you think.

[neich]

Hi @zhouzx-UCAS , sorry for the delay.

I see your point about clipping. Maybe we should implement a mechanism so each model has a kind of validate_state() that after computing a new state, it makes all clipping and any other check necessary, so all the next computations (including coupling) use a valid state. Would that be acceptable?

About, the FIC issue, yes, your contribution is definitely welcome. Make a separate PR and I will accept it

[zhouzx-UCAS]

Hi @neich,

Thank you very much for the reply.

I will make a separate PR for the FIC issue.

And yes, I think a validate_state() mechanism would be a good solution here. If you are okay with it, I can try implementing it and submit a PR for your review.

Best,
Zi-Xuan

[neich]

@zhouzx-UCAS I have created PR #109 to deal with the generic validation of state variables. Check it and if it looks reasonable I will merge it

[neich]

@zhouzx-UCAS I think that the best course of action is just close this PR, then approve #109 and then you can open a new one with the FIC fixes

[zhouzx-UCAS]

Hi @neich,
I'm testing PR #109 and will let you know the results soon!

[zhouzx-UCAS]

Hi @neich,

I have now carefully read through your PR #109, and I think the logic is sound.

In my tests, the new implementation reproduces the trend reported in the paper, and I found that there are small numerical differences (under the same random seed) compared with the results I previously obtained by clipping directly inside LinearCouplingModel.
I am not 100% sure about the exact source of that difference, but my current impression is that it may come from the fact that PR #109 projects the state itself back into the valid range (which is the more rigorous approach and obviously correct), whereas in my earlier modification I only clipped the value used in coupling, without writing the clipped state back. As a result, the next update step will still be using the unclipped state.

I have also updated the Deco2018 model to be compatible with the new PR #109 mechanism and submitted that as PR #110. In addition, I submitted the FIC fix separately as PR #111. Please let me know if you have any concerns or thoughts about these PRs.

Best,
Zi-Xuan