Updated coupling and other minor things!

Projects/AD-ABeta_and_Tau/DMF_AD.py

@@ -0,0 +1,179 @@
+# ==========================================================================
+# ==========================================================================
+# ==========================================================================
+# Model for AD simulations.
+#
+# Described at:
+# Whole-brain modeling of the differential influences of Amyloid-Beta and Tau in Alzheimer`s Disease, by
+# Gustavo Patow, Leon Stefanovski, Petra Ritter, Gustavo Deco, Xenia Kobeleva and for the
+# Alzheimer’s Disease Neuroimaging Initiative. Accepted at Alzheimer's Research & Therapy, 2023
+#
+# Adapted Dynamic Mean Field (DMF) model (a.k.a., Reduced Wong-Wang), from
+#
+# [Deco_2014] G. Deco, A. Ponce-Alvarez, P. Hagmann, G.L. Romani, D. Mantini, M. Corbetta
+#             How local excitation-inhibition ratio impacts the whole brain dynamics
+#             J. Neurosci., 34 (2014), pp. 7886-7898
+#
+# Modifyied by Gustavo Patow for Alzheimer's Disease simulations
+# ==========================================================================
+import numpy as np
+from numba import jit
+
+print("Going to use the AD Dynamic Mean Field (adDMF) neuronal model...")
+
+
+def recompileSignatures():
+    # Recompile all existing signatures. Since compiling isn’t cheap, handle with care...
+    # However, this is "infinitely" cheaper than all the other computations we make around here ;-)
+    phie.recompile()
+    phii.recompile()
+    dfun.recompile()
+
+
+# ==========================================================================
+# ==========================================================================
+# ==========================================================================
+# Dynamic Mean Field (DMF) Model Constants
+# --------------------------------------------------------------------------
+taon = 100.
+taog = 10.
+gamma_e = 0.641 / 1000
+gamma_i = 1. / 1000.
+J_NMDA = 0.15       # [nA] NMDA current
+I0 = 0.382  #.397  # [nA] overall effective external input
+Jexte = 1.
+Jexti = 0.7
+w = 1.4
+we = 2.1        # Global coupling scaling (G in the paper)
+SC = None       # Structural connectivity (should be provided externally)
+
+
+# ==========================================================================
+# AD param initialization...
+# --------------------------------------------------------------------------
+M_e = 1.    # WARNING: In general, this must be initialized outside!
+M_i = 1.
+Abeta = None
+Tau = None
+
+def set_AD_Burden(coefs):
+    global M_e, M_i
+    # First, the excitatory contribution
+    M_e = (1 + coefs[0] + coefs[1] * Abeta) * (1 + coefs[2] + coefs[3] * Tau)
+    # Now, the inhibitory contribution
+    M_i = (1 + coefs[4] + coefs[5] * Abeta) * (1 + coefs[6] + coefs[7] * Tau)
+
+    recompileSignatures()
+
+
+# --------------------------------------------------------------------------
+# Simulation variables
+# @jit(nopython=True)
+def initSim(N):
+    sn_sg = 0.001 * np.ones((2, N))  # Here we initialize with np.ones, but other models use np.zeros...
+    return sn_sg
+
+
+J = None    # WARNING: In general, J must be initialized outside!
+def initJ(N):  # A bit silly, I know...
+    global J
+    J = np.ones(N)
+
+
+# --------------------------------------------------------------------------
+# Variables of interest, needed for bookkeeping tasks...
+# xn = None  # xn, excitatory synaptic activity
+# rn = None  # rn, excitatory firing rate
+# @jit(nopython=True)
+def numObsVars():  # Returns the number of observation vars used, here xn and rn
+    return 3
+
+
+# --------------------------------------------------------------------------
+# Set the parameters for this model
+def setParms(modelParms):
+    global we, J, SC, M_e, M_i
+    if 'we' in modelParms:
+        we = modelParms['we']
+    if 'J' in modelParms:
+        J = modelParms['J']
+    if 'SC' in modelParms:
+        SC = modelParms['SC']
+    if 'M_e' in modelParms:
+        M_e = modelParms['M_e']
+    if 'M_i' in modelParms:
+        M_i = modelParms['M_i']
+    if 'coefs' in modelParms:
+        set_AD_Burden(modelParms('coefs'))
+
+
+def getParm(parmList):
+    if 'we' in parmList:
+        return we
+    if 'J' in parmList:
+        return J
+    if 'be' in parmList:
+        return be
+    if 'ae' in parmList:
+        return ae
+    if 'SC' in parmList:
+        return SC
+    return None
+
+
+# -----------------------------------------------------------------------------
+# ----------------- Dynamic Mean Field (a.k.a., reducedWongWang) --------------
+# -----------------------------------------------------------------------------
+
+# ----------------- Coumpling ----------------------
+# @jit(nopython=True)
+# def instantaneousDirectCouplingOp(x):
+#     return SC @ x      # coupling = SC @ sn
+couplingOp = None
+
+
+# --------------------------------------------------------------------------
+# transfer WholeBrain:
+# --------------------------------------------------------------------------
+# transfer function: excitatory
+ae=310.
+be=125.
+de=0.16
+@jit(nopython=True)
+def phie(x):
+    y = M_e * (ae*x-be)  # M_e will be set by function set_AD_Burden
+    return y/(1.-np.exp(-de*y))
+
+
+# transfer function: inhibitory
+ai=615.
+bi=177.
+di=0.087
+@jit(nopython=True)
+def phii(x):
+    y = M_i * (ai*x-bi)  # This is for Modality A in our study...
+                         # M_i will be set by function set_AD_Burden
+    return y/(1.-np.exp(-di*y))
+
+
+# transfer WholeBrain used by the simulation...
+He = phie
+Hi = phii
+
+
+@jit(nopython=True)
+def dfun(simVars, coupling, I_external):
+    sn = simVars[0]; sg = simVars[1]  # should be [sn, sg] = simVars
+    # This is the regular DMF behaviour
+    coupl = coupling.couple(sn)
+    xn = I0 * Jexte + w * J_NMDA * sn + we * J_NMDA * coupl - J * sg + I_external  # Eq for I^E (5). I_external = 0 => resting state condition.
+    xg = I0 * Jexti + J_NMDA * sn - sg  # Eq for I^I (6). \lambda = 0 => no long-range feedforward inhibition (FFI)
+    rn = He(xn)  # Calls He(xn). r^E = H^E(I^E) in the paper (7)
+    rg = Hi(xg)  # Calls Hi(xg). r^I = H^I(I^I) in the paper (8)
+    dsn = -sn / taon + (1. - sn) * gamma_e * rn
+    dsg = -sg / taog + rg * gamma_i
+    return np.stack((dsn, dsg)), np.stack((xn, rn, rg))
+
+# ==========================================================================
+# ==========================================================================
+# ==========================================================================EOF

Projects/AD-ABeta_and_Tau/Prepro.py

@@ -0,0 +1,135 @@
+# --------------------------------------------------------------------------------------
+# Full pipeline for processing AD, MCI and HC subjects
+#
+# Described at:
+# Whole-brain modeling of the differential influences of Amyloid-Beta and Tau in Alzheimer`s Disease, by
+# Gustavo Patow, Leon Stefanovski, Petra Ritter, Gustavo Deco, Xenia Kobeleva and for the
+# Alzheimer’s Disease Neuroimaging Initiative. Accepted at Alzheimer's Research & Therapy, 2023
+#
+# By Gustavo Patow
+#
+# Pre-requisites:
+#   Before executing this, be sure to have correctly configured the setup_AD.py file...
+#
+# --------------------------------------------------------------------------------------
+import numpy as np
+import scipy.io as sio
+
+from setup import *
+
+import WholeBrain.Observables.FC as FC
+import WholeBrain.Observables.swFCD as swFCD
+import WholeBrain.Observables.phFCD as phFCD
+from WholeBrain.Utils.plotSC import plotSC
+
+
+# =====================================================================
+# =====================================================================
+#                      Single Subject Pipeline
+# =====================================================================
+# =====================================================================
+def preprocessingPipeline(subject, SCnorm, all_fMRI,  #, abeta,
+                          distanceSettings,  # This is a dictionary of {name: (distance module, apply filters bool)}
+                          WEs,  # wStart=0.0, wEnd=6.0, wStep=0.05,
+                          plotMaxFrecForAllWe=False):
+    fileName = save_folder + '/FICWeights-'+subject+'/BenjiBalancedWeights-{}.mat'
+
+    print("###################################################################\n"*2+"#")
+    print(f"## Pre-processing pipeline on {subject}!!!\n#")
+    print("###################################################################\n"*2+"\n")
+    print(f"# Compute BalanceFIC ({wStart} to {wEnd} with step {wStep})")
+    print("###################################################################")
+    # BalanceFIC.useDeterministicIntegrator = useDeterministicIntegrator
+    BalanceFIC.verbose = True
+    balancedParms = BalanceFIC.Balance_AllJ9(SCnorm, WEs=WEs,  # wStart=wStart, wEnd=wEnd, wStep=wStep,
+                                             baseName=fileName)
+    modelParms = [balancedParms[i] for i in balancedParms]
+
+    # Let's plot it as a verification measure...
+    if plotMaxFrecForAllWe:
+        import Tests.DecoEtAl2014.fig2 as Fig2
+        Fig2.plotMaxFrecForAllWe(SCnorm, wStart=WEs[0], wEnd=WEs[-1], wStep=WEs[1]-WEs[0],
+                                 extraTitle='', precompute=False, fileName=fileName)  # We already precomputed everything, right?
+
+    # Now, optimize all we (G) values: determine optimal G to work with
+    print("\n\n###################################################################")
+    print("# Compute ParmSweep")
+    print("###################################################################\n")
+    outFilePath = save_folder+'/'+subject+'-temp'
+    fitting = ParmSeep.distanceForAll_Parms(all_fMRI, WEs, modelParms, NumSimSubjects=len(all_fMRI),
+                                            distanceSettings=distanceSettings,
+                                            parmLabel='we',
+                                            outFilePath=outFilePath)
+
+    optimal = {sd: distanceSettings[sd][0].findMinMax(fitting[sd]) for sd in distanceSettings}
+    return optimal
+
+
+# =====================================================================
+# =====================================================================
+#                            main
+# =====================================================================
+# =====================================================================
+def processRangeValues(argv):
+    import getopt
+    try:
+        opts, args = getopt.getopt(argv,'',["wStart=","wEnd=","wStep="])
+    except getopt.GetoptError:
+        print('Prepro.py --wStart <wStartValue> --wEnd <wEndValue> --wStep <wStepValue>')
+        sys.exit(2)
+    wStart = 0.; wEnd = 5.5; wStep = 0.05
+    for opt, arg in opts:
+        if opt == '-h':
+            print('Prepro.py -wStart <wStartValue> -wEnd <wEndValue> -wStep <wStepValue>')
+            sys.exit()
+        elif opt in ("--wStart"):
+            wStart = float(arg)
+        elif opt in ("--wEnd"):
+            wEnd = float(arg)
+        elif opt in ("--wStep"):
+            wStep = float(arg)
+    print(f'Input values are: wStart={wStart}, wEnd={wEnd}, wStep={wStep}')
+    return wStart, wEnd, wStep
+
+
+visualizeAll = True
+if __name__ == '__main__':
+    import sys
+    wStart, wEnd, wStep = processRangeValues(sys.argv[1:])  # Default: wStart = 0.; wEnd = 5.5; wStep = 0.05
+
+    # --------------------------------------------------
+    # Compute the average SC for the HC subjects
+    # --------------------------------------------------
+    avgSCMatrix = dataLoader.computeAvgSC_HC_Matrix(classification, base_folder + "/connectomes")
+    dataLoader.analyzeMatrix("AvgHC", avgSCMatrix)
+    finalAvgMatrixHC = dataLoader.correctSC(avgSCMatrix)
+    sio.savemat(save_folder + '/AvgHC_SC.mat', {'SC':finalAvgMatrixHC})
+    dataLoader.analyzeMatrix("AvgHC norm", finalAvgMatrixHC)
+    print("# of elements in AVG connectome: {}".format(finalAvgMatrixHC.shape))
+    # plotSC.justPlotSC('AVG<HC>', finalMatrix, plotSC.plotSCHistogram)
+    # plot_cc_empSC_empFC(HCSubjects)
+    neuronalModel.setParms({'SC': finalAvgMatrixHC})
+    neuronalModel.couplingOp = Couplings.instantaneousDirectCoupling(finalAvgMatrixHC)
+
+    # --------------------------------------------------
+    # load all HC fMRI data
+    # --------------------------------------------------
+    all_HC_fMRI = dataLoader.load_fullCohort_fMRI(classification, base_folder)
+
+    # Configure and compute Simulation
+    # ------------------------------------------------
+    subjectName = 'AvgHC'
+    distanceSettings = {'FC': (FC, False), 'swFCD': (swFCD, True), 'phFCD': (phFCD, True)}
+
+    BalanceFIC.use_N_algorithm = False  # To make sure we use Gus' algorithm
+    WEs = np.arange(wStart, wEnd+0.01, wStep)
+    optimal = preprocessingPipeline(subjectName, finalAvgMatrixHC, all_HC_fMRI,
+                                    distanceSettings,
+                                    WEs=WEs,  # wStart=wStart, wEnd=wEnd+0.01, wStep=wStep,
+                                    plotMaxFrecForAllWe=False)
+    print (f"Last info: Optimal in the CONSIDERED INTERVAL only: {wStart}, {wEnd}, {wStep} (not in the whole set of results!!!)")
+    print("".join(f" - Optimal {k}({optimal[k][1]})={optimal[k][0]}\n" for k in optimal))
+
+# ================================================================================================================
+# ================================================================================================================
+# ================================================================================================================EOF

Projects/AD-ABeta_and_Tau/README.md

@@ -1,13 +1,21 @@
 # Whole-brain modeling of the differential influences of Amyloid-Beta and Tau in Alzheimer`s Disease
 
-Python code for the paper
+Python code for the paper:
 
 Whole-brain modeling of the differential influences of Amyloid-Beta and Tau in Alzheimer`s Disease, by 
 Gustavo Patow, Leon Stefanovski, Petra Ritter, Gustavo Deco, Xenia Kobeleva and for the
 Alzheimer’s Disease Neuroimaging Initiative. Accepted at Alzheimer's Research & Therapy, 2023
 
-This code is "broken research code", and thus the following disclaimer appllies:
-THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
-WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
-COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
-OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
+This code is "broken research code", and thus the following disclaimer appllies: 
+THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+## Code usage
+1) setup.py is the main configuration file. It should be modified first. This file uses dataLoader.py to perform the actual loading of the files.
+2) Prepro.py should be executed next, as it computes the preprocessing of the model (computation of G, we in the code) and uses the previous one.
+3) plotPrepro.py plots the results of Prepro.py. Should be done AFTER Prepro.py.
+4) Finally, pipeline2.py should be executed. This file was written before the big re-factorization of the WholeBrain library, so it partially uses its newer features. Also, it uses the G=3.1 value found at Prepro.py.
+5) plot3DBrainBurdens is used to plot the nice 3D plots ocmparing ABeta and Tau for HC, MCI and AD.
+6) Shuffling.py is used for the comparisons with shuffled random data...
+7) plotParmComparisonAcrossCohort.py is used for plotting the different figures in the paper, see comments within at the respective boolean variables...
+
+If the AT(N) classification should be used, first run classific_AT(N).py, and then configure and use setup_ATN.py. All the rest of the system should continue working as before (this only affects the final cohort grouping for plotting).