OSCR

Sensitivity analysis enlightens effects of connectivity in a Neural Mass Model under Control-Target mode.

Code ↔ Paper

The paper beside its authors' code: matches between them have not been computed for this paper yet.

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 780 lines · 19 KB · no license

  1. #!/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3. """
  4. Created on Tue Jan 20 14:39:12 2026
  5. @author: anais
  6. """
  7. from numpy import *
  8. import math
  9. #%% parameters for the neural model
  10. def get_FixedParametersForCurrents() :
  11. Jnmda=0.15
  12. G=0.69
  13. Wplus = 1.4
  14. Wminus = 1
  15. Wi = 0.7
  16. We = 1
  17. I0 = 0.382
  18. global autocoupling_E
  19. global basal_forcing_E
  20. global coupling_EtoI
  21. global autocoupling_I_abs
  22. global basal_forcing_I
  23. autocoupling_E = Wplus*Jnmda
  24. basal_forcing_E = We * I0
  25. coupling_EtoI = Jnmda
  26. autocoupling_I_abs = Wminus #inhibitory autocoupling (in absolute value)
  27. basal_forcing_I = Wi*I0
  28. return(autocoupling_E,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,Jnmda,G)
  29. def get_FixedParametersForStates() :
  30. BetaE=6.6
  31. AlphaE=72
  32. BetaI=180
  33. AlphaI=530
  34. return(BetaE,AlphaE,BetaI,AlphaI)
  35. def get_FixedParametersForExcFiringRates() :
  36. a=310.
  37. b=125
  38. d=0.16
  39. return(a,b,d)
  40. def get_FixedParametersForInhFiringRates() :
  41. a=615.
  42. b=177
  43. d=0.087
  44. return(a,b,d)
  45. #%% parameters for the hemodynamic and bold models
  46. def get_FixedParametersForHemodynamics() :
  47. alpha=0.32
  48. gamma=0.41
  49. kappa=0.65
  50. tau=0.98
  51. rho=0.34
  52. return(alpha,gamma,kappa,tau,rho)
  53. def get_FixedParametersForBold() :
  54. Theta0=28.265*3
  55. rho=0.34
  56. TE=0.03 #30ms
  57. epsilon=0.47
  58. r0=110
  59. V0=0.02
  60. k1 = 4.3*Theta0*rho*TE
  61. k2 = epsilon*r0*rho*TE
  62. k3 = 1-epsilon
  63. return(V0,k1,k2,k3)
  64. #%% Fixing the type of connectivity
  65. def get_Matrix_kE_ForATypeOfConnectivity(TypeOfConnectivity):
  66. if TypeOfConnectivity=='EE':
  67. k = array([[0., 1.],
  68. [1., 0.]])
  69. elif TypeOfConnectivity=='II':
  70. k = array([[0., 0.],
  71. [0., 0.]])
  72. elif TypeOfConnectivity=='EI':
  73. k = array([[0., 1.],
  74. [0., 0.]])
  75. elif TypeOfConnectivity == 'IE':
  76. k = array([[0., 0.],
  77. [1., 0.]])
  78. return(k)
  79. def get_k_inh(k): # for the inhibitory pool. k_inh is a matrix that contains 0 on the diagonal and (1-k) everywhere else.
  80. mat1_0diago = ones(k.shape)
  81. for i in range(len(k)):
  82. mat1_0diago[i,i] = 0 #mat1_0diago contains 0 on the diagonal and 1 everywhere else
  83. k_inh = mat1_0diago-k
  84. return(k_inh)
  85. #%% Simulations with the neural model
  86. def FindFixedPoint(autocoupling_E,coupling_ItoE_abs,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,M,G,k,Jnmda,Tglu,Tgaba,N): #cf fonction EtudeStabilite papier tuning
  87. #----------Duration of the first simulation and timeline
  88. T0=0
  89. Tmax = 10
  90. (Time,dt) = get_Timeline(T0,Tmax)
  91. #---- Stimulation forcing on the excitatory pools
  92. StimulationFunction = zeros((N,len(Time)))
  93. #------- Initial conditions
  94. dsn_dt = 1
  95. dsg_dt = 1
  96. sn_fp = 0.5*ones((N,1))
  97. sg_fp = 0.5*ones((N,1))
  98. sn_fp = 0.1*ones((N,1))
  99. sg_fp = 0.05*ones((N,1))
  100. simu=0
  101. while all(abs(dsn_dt)>1e-10) and all(abs(dsg_dt)>1e-10): #We continue the simulation until we get close enough to the fixed point
  102. simu=simu+1
  103. Sn0 = sn_fp
  104. Sg0 = sg_fp
  105. (Xn,Xg,Rn,Rg,dSndt,dSgdt,Sn,Sg) = get_Dynamics(autocoupling_E,coupling_ItoE_abs,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,M,G,k,Jnmda,Tglu,Tgaba,Sn0,Sg0,dt,Time,StimulationFunction)
  106. sn_fp = reshape(Sn[:,-1],(N,1))
  107. sg_fp = reshape(Sg[:,-1],(N,1))
  108. dsn_dt = reshape(dSndt[:,-1],(N,1))
  109. dsg_dt = reshape(dSgdt[:,-1],(N,1))
  110. xn_fp = reshape(Xn[:,-1],(N,1))
  111. xg_fp = reshape(Xg[:,-1],(N,1))
  112. rn_fp = reshape(Rn[:,-1],(N,1))
  113. rg_fp = reshape(Rg[:,-1],(N,1))
  114. return(sn_fp,sg_fp,xn_fp,xg_fp,rn_fp,rg_fp)
  115. def get_Timeline(T0,Tmax):
  116. dt=0.001
  117. sampling=round((Tmax-T0)/dt)
  118. Tmax = T0 + dt*sampling
  119. Time=linspace(T0,Tmax,int(sampling)+1)
  120. return(Time,dt)
  121. def get_Dynamics(autocoupling_E,coupling_ItoE_abs,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,M,G,k,Jnmda,Tglu,Tgaba,Sn0,Sg0,dt,Time,StimulationFunction):
  122. N= Sn0.shape[0]
  123. # memory allocation
  124. Sn = zeros((N,len(Time)))
  125. Sg = zeros((N,len(Time)))
  126. dSndt = zeros((N,len(Time)))
  127. dSgdt = zeros((N,len(Time)))
  128. Xn = zeros((N,len(Time)))
  129. Xg = zeros((N,len(Time)))
  130. Rn = zeros((N,len(Time)))
  131. Rg = zeros((N,len(Time)))
  132. # we save the initial condition
  133. Sn[:,0] = reshape(Sn0,(Sn0.T).shape)
  134. Sg[:,0] = reshape(Sg0,(Sg0.T).shape)
  135. for j in range(len(Time)):
  136. sn_past = reshape(Sn[:,j],(N,1))
  137. sg_past = reshape(Sg[:,j],(N,1))
  138. t=Time[j]
  139. (xn,xg,rn,rg,dsndt,dsgdt,sn_new,sg_new) = Integration(autocoupling_E,coupling_ItoE_abs,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,M,G,k,Jnmda,Tglu,Tgaba,sn_past,sg_past,dt,j,StimulationFunction)
  140. dSndt[:,j] = reshape(dsndt,(dsndt.T).shape)
  141. dSgdt[:,j] = reshape(dsgdt,(dsgdt.T).shape)
  142. Rn[:,j] = reshape(rn,(rn.T).shape)
  143. Rg[:,j] = reshape(rg,(rg.T).shape)
  144. Xn[:,j] = reshape(xn,(xn.T).shape)
  145. Xg[:,j] = reshape(xg,(xg.T).shape)
  146. if j<len(Time)-1: # We do not keep the last point (i.e. Sn(t+dt) obtained with dSndt(t), for t=Time[-1]). Thus, Sn and dSndt have the same dimensions and the last point of each array corresponds to Sn(t) and dSndt(t) with t=Time[-1]
  147. Sn[:,j+1] = reshape(sn_new,(sn_new.T).shape)
  148. Sg[:,j+1] = reshape(sg_new,(sg_new.T).shape)
  149. return(Xn,Xg,Rn,Rg,dSndt,dSgdt,Sn,Sg)
  150. def Integration(autocoupling_E,coupling_ItoE_abs,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,M,G,k,Jnmda,Tglu,Tgaba,sn,sg,dt,j,StimulationFunction):
  151. StimulationFunction_at_t = reshape(StimulationFunction[:,j],sn.shape)
  152. (xn,xg,rn,rg,dsndt,dsgdt) = get_Derivatives(sn,sg,autocoupling_E,coupling_ItoE_abs,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,M,G,k,Jnmda,Tglu,Tgaba,j,StimulationFunction_at_t)
  153. N = sn.shape[0]
  154. sn_next=zeros(sn.shape, dtype=float)
  155. sg_next=zeros(sg.shape, dtype=float)
  156. sn_next = sn + dt*dsndt
  157. sg_next = sg + dt*dsgdt
  158. return(xn,xg,rn,rg,dsndt,dsgdt,sn_next,sg_next)
  159. def get_Derivatives(sn,sg,autocoupling_E,coupling_ItoE_abs,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,M,G,k,Jnmda,Tglu,Tgaba,j,StimulationFunction_at_t):
  160. #--- Currents
  161. xn = ExcitatoryCurrent(autocoupling_E,coupling_ItoE_abs,basal_forcing_E,sn,sg,j,StimulationFunction_at_t,M,G,k,Jnmda)
  162. xg = InhibitoryCurrent (coupling_EtoI,autocoupling_I_abs,basal_forcing_I,sn,sg,j,M,G,k,Jnmda)
  163. #--- Firing rates
  164. rn=ExcitatoryFiringRate(xn)
  165. rg=InhibitoryFiringRate(xg)
  166. #--- Differential equations on state variables
  167. (BetaE,AlphaE,BetaI,AlphaI) = get_FixedParametersForStates()
  168. dsndt = -BetaE*sn+AlphaE*Tglu*(1-sn)*rn
  169. dsgdt = -BetaI*sg+AlphaI*Tgaba*(1-sg)*rg
  170. return(xn,xg,rn,rg,dsndt,dsgdt)
  171. #--- functions for firing rates
  172. def ExcitatoryFiringRate(x):
  173. (a,b,d) = get_FixedParametersForExcFiringRates()
  174. y=a*x-b
  175. result=zeros(y.shape, dtype=float)
  176. for i in range(len(y)):
  177. result[i] = y[i]/(1-math.exp(-d*y[i]))
  178. # if y[i]!=0:
  179. # result[i] = y[i]/(1-math.exp(-d*y[i]))
  180. # else:
  181. # result[i]=0
  182. return(result)
  183. def InhibitoryFiringRate(x):
  184. (a,b,d) = get_FixedParametersForInhFiringRates()
  185. y=a*x-b
  186. result=zeros(y.shape, dtype=float)
  187. for i in range(len(y)):
  188. result[i] = y[i]/(1-math.exp(-d*y[i]))
  189. # if y[i]!=0:
  190. # result[i] = y[i]/(1-math.exp(-d*y[i]))
  191. # else:
  192. # result[i]=0
  193. return(result)
  194. #--- functions for currents
  195. def ExcitatoryCurrent (autocoupling_E,coupling_ItoE_abs,basal_forcing_E,sn,sg,j,StimulationFunction_at_t,M,G,k,Jnmda):
  196. xn = autocoupling_E*sn - coupling_ItoE_abs*sg + basal_forcing_E + StimulationFunction_at_t + G*Jnmda*dot( k*M,sn )
  197. return(xn)
  198. def InhibitoryCurrent (coupling_EtoI,autocoupling_I_abs,basal_forcing_I,sn,sg,j,M,G,k,Jnmda):
  199. k_inh = get_k_inh(k)
  200. xg = coupling_EtoI*sn - autocoupling_I_abs*sg + basal_forcing_I + G*Jnmda*dot( k_inh*M,sn )
  201. return(xg)
  202. #%% Simulations with the hemodynamic and bold models
  203. def HemoBoldDynamics(SnDyna,z0,f0,v0,q0,dt,Time):
  204. bold0=BoldFunction(v0,q0)
  205. N= SnDyna.shape[0]
  206. #allocation mémoire
  207. Z = zeros((N,len(Time)))
  208. F = zeros((N,len(Time)))
  209. V = zeros((N,len(Time)))
  210. Q = zeros((N,len(Time)))
  211. BOLD = zeros((N,len(Time)))
  212. # we keep the initial condition
  213. Z[:,0] = reshape(z0,(z0.T).shape)
  214. F[:,0] = reshape(f0,(f0.T).shape)
  215. V[:,0] = reshape(v0,(v0.T).shape)
  216. Q[:,0] = reshape(q0,(q0.T).shape)
  217. BOLD[:,0] = reshape(bold0,(bold0.T).shape)
  218. for j in range(1,len(Time)):
  219. sn_past = SnDyna[0:N,j-1:j]
  220. z_past = reshape(Z[:,j-1],(N,1))
  221. f_past = reshape(F[:,j-1],(N,1))
  222. v_past = reshape(V[:,j-1],(N,1))
  223. q_past = reshape(Q[:,j-1],(N,1))
  224. bold_past=BoldFunction(v_past,q_past)
  225. (z_new,f_new,v_new,q_new) = HemodynamicsIntegration(sn_past,z_past,f_past,v_past,q_past,dt)
  226. bold_new=BoldFunction(v_new,q_new)
  227. Z[:,j] = reshape(z_new,(z_new.T).shape)
  228. F[:,j] = reshape(f_new,(f_new.T).shape)
  229. V[:,j] = reshape(v_new,(v_new.T).shape)
  230. Q[:,j] = reshape(q_new,(q_new.T).shape)
  231. BOLD[:,j] = reshape(bold_new,(bold_new.T).shape)
  232. return(Z,F,V,Q,BOLD)
  233. def HemodynamicsIntegration(sn,z,f,v,q,dt):
  234. (alpha,gamma,kappa,tau,rho) = get_FixedParametersForHemodynamics()
  235. dzdt = sn-kappa*z-gamma*(f-1)
  236. dfdt = z
  237. dvdt = (1/tau)*(f-pow(v,1/alpha))
  238. dqdt = (1/(tau*rho))*f*(1-pow((1-rho),1/f)) - (1/tau)*q*pow(v,(1/alpha)-1)
  239. z_next=zeros(z.shape, dtype=float)
  240. f_next=zeros(f.shape, dtype=float)
  241. v_next=zeros(v.shape, dtype=float)
  242. q_next=zeros(q.shape, dtype=float)
  243. z_next=z+dt*dzdt
  244. f_next=f+dt*dfdt
  245. v_next=v+dt*dvdt
  246. q_next=q+dt*dqdt
  247. return(z_next,f_next,v_next,q_next)
  248. def BoldFunction(v,q):
  249. (V0,k1,k2,k3) = get_FixedParametersForBold()
  250. y = V0*(k1*(1-q)+k2*(1-q*(1/v))+k3*(1-v))
  251. return(y)
  252. #%%% Sensitivities
  253. #%% Sensitivities of the isolated pools (with their autocoupling)
  254. def CloseLoopSensitivity_Pool_E(Sn,zE,autocoupling_E,Tglu):
  255. Sensi_Sn_zE = zeros(Sn.shape)
  256. (BetaE,AlphaE,BetaI,AlphaI) = get_FixedParametersForStates()
  257. (XN,RN,DRN_DXN) = Calculate_XN_RN_DRNDXN_PoolE(Sn,autocoupling_E,zE)
  258. # Sensitivity of sn w.r.t. the forcing zE
  259. Sensi_Sn_zE = 1/( ( (BetaE + AlphaE*Tglu*RN)/(AlphaE*Tglu*(1-Sn)*DRN_DXN) ) - autocoupling_E)
  260. # Feedback gain
  261. Gain = (AlphaE*Tglu*(1-Sn)*DRN_DXN*autocoupling_E)/(BetaE + AlphaE*Tglu*RN)
  262. return(Sensi_Sn_zE, Gain)
  263. def CloseLoopSensitivity_Pool_I(Sg,zI,autocoupling_I_abs,Tgaba):
  264. Sensi_Sg_zI = zeros(Sg.shape)
  265. (BetaE,AlphaE,BetaI,AlphaI) = get_FixedParametersForStates()
  266. (XG,RG,DRG_DXG) = Calculate_XG_RG_DRGDXG_PoolI(Sg,autocoupling_I_abs,zI)
  267. # Sensitivity of sg w.r.t. the forcing zI
  268. Sensi_Sg_zI = 1/( ( (BetaI + AlphaI*Tgaba*RG)/(AlphaI*Tgaba*(1-Sg)*DRG_DXG) ) + autocoupling_I_abs )
  269. # Feedback gain
  270. Gain = (AlphaI*Tgaba*(1-Sg)*DRG_DXG * (-autocoupling_I_abs))/(BetaI + AlphaI*Tgaba*RG)
  271. return(Sensi_Sg_zI, Gain)
  272. def Calculate_XN_RN_DRNDXN_PoolE(Sn,autocoupling_E,zE):
  273. XN = zeros(Sn.shape)
  274. RN = zeros(Sn.shape)
  275. DRN_DXN = zeros(Sn.shape)
  276. (Aexc,Bexc,Dexc) = get_FixedParametersForExcFiringRates()
  277. XN = autocoupling_E*Sn + zE
  278. RN = (Aexc * XN-Bexc) /(1-exp(-Dexc*(Aexc * XN - Bexc)))
  279. DRN_DXN = (Aexc/2) * ( (exp(Dexc*(Aexc * XN - Bexc)) -1 - Dexc*(Aexc * XN - Bexc) ) / (cosh(Dexc*(Aexc * XN - Bexc)) -1) )
  280. return(XN,RN,DRN_DXN)
  281. def Calculate_XG_RG_DRGDXG_PoolI(Sg,autocoupling_I_abs,zI):
  282. XG = zeros(Sg.shape)
  283. RG = zeros(Sg.shape)
  284. DRG_DXG = zeros(Sg.shape)
  285. (Ainh,Binh,Dinh) = get_FixedParametersForInhFiringRates()
  286. XG = - autocoupling_I_abs*Sg + zI
  287. RG = (Ainh * XG - Binh) /(1-exp(-Dinh*(Ainh * XG - Binh)))
  288. DRG_DXG =(Ainh/2) * ( (exp(Dinh*(Ainh * XG - Binh)) -1 - Dinh*(Ainh * XG - Binh) ) / (cosh(Dinh*(Ainh * XG - Binh)) -1) )
  289. return(XG,RG,DRG_DXG)
  290. #%% Sensitivities of the isolated regions (as functions of open loop sensitivities of the pools)
  291. #---- Open loop sensitivities of the pools
  292. def OpenLoopSensitivities_Pools(Sn,Sg,zE,zI,autocoupling_E,autocoupling_I_abs,Tglu,Tgaba):
  293. (Sensi_Sn_zE, Gain_E) = CloseLoopSensitivity_Pool_E(Sn,zE,autocoupling_E,Tglu)
  294. (Sensi_Sg_zI, Gain_I) = CloseLoopSensitivity_Pool_I(Sg,zI,autocoupling_I_abs,Tgaba)
  295. return(Sensi_Sn_zE, Sensi_Sg_zI)
  296. #---- Close loop sensitivities of the isolated area
  297. def CloseLoopSensitivities_OneRegion(Sn,Sg,xE,xI,autocoupling_E,coupling_ItoE_abs,coupling_EtoI,autocoupling_I_abs,Tglu,Tgaba):
  298. (BetaE,AlphaE,BetaI,AlphaI) = get_FixedParametersForStates()
  299. # ---initialisation
  300. # open loop sensitivities
  301. Sensi_Sn_xE_O = zeros(Sn.shape)
  302. Sensi_Sg_xI_O = zeros(Sn.shape)
  303. # close loop sensitivities
  304. Sensi_Sn_xE = zeros(Sn.shape) # sensitivity of sn w.r.t. the forcing xE
  305. Sensi_Sg_xE = zeros(Sn.shape) # sensitivity of sg w.r.t. the forcing xE
  306. Sensi_Sn_xI = zeros(Sn.shape) # sensitivity of sn w.r.t. the forcing xI
  307. Sensi_Sg_xI = zeros(Sn.shape) # sensitivity of sg w.r.t. the forcing xI
  308. # We evaluate the open loop sensitivities using the fixed point in the close loop system
  309. ZE = - coupling_ItoE_abs*Sg + xE
  310. ZI = coupling_EtoI*Sn + xI
  311. (Sensi_Sn_xE_O, Sensi_Sg_xI_O) = OpenLoopSensitivities_Pools(Sn,Sg,ZE,ZI,autocoupling_E,autocoupling_I_abs,Tglu,Tgaba)
  312. # Feedback gain
  313. Gain = -coupling_EtoI * coupling_ItoE_abs * Sensi_Sn_xE_O * Sensi_Sg_xI_O
  314. # We calculate the close loop sensitivities, using the open loop sensitivities and the feedback gain
  315. Sensi_Sn_xE = Sensi_Sn_xE_O/(1-Gain)
  316. Sensi_Sg_xE = (coupling_EtoI * Sensi_Sn_xE_O * Sensi_Sg_xI_O )/(1-Gain)
  317. Sensi_Sn_xI = (-coupling_ItoE_abs * Sensi_Sn_xE_O * Sensi_Sg_xI_O)/(1-Gain)
  318. Sensi_Sg_xI = Sensi_Sg_xI_O/(1-Gain)
  319. return(Sensi_Sn_xE,Sensi_Sg_xE,Sensi_Sn_xI,Sensi_Sg_xI,Sensi_Sn_xE_O,Sensi_Sg_xI_O,Gain)
  320. #%% Sensitivities of two regions (as functions of open loop sensitivitie of each region)
  321. #---- Open loop sensitivities of one region
  322. def OpenLoopSensitivities_OneRegion(Sni,Sgi,xEi,xIi,autocoupling_E_i,coupling_ItoE_abs_i,coupling_EtoI_i,autocoupling_I_abs_i,Tglui,Tgabai): #give the input arguments only for the region of interest
  323. (Sensi_Sni_xEi,Sensi_Sgi_xEi,Sensi_Sni_xIi,Sensi_Sgi_xIi,Sensi_Sni_xEi_O,Sensi_Sgi_xIi_O,Gain) = CloseLoopSensitivities_OneRegion(Sni,Sgi,xEi,xIi,autocoupling_E_i,coupling_ItoE_abs_i,coupling_EtoI_i,autocoupling_I_abs_i,Tglui,Tgabai)
  324. return(Sensi_Sni_xEi,Sensi_Sgi_xEi,Sensi_Sni_xIi,Sensi_Sgi_xIi)
  325. #---- Close loop sensitivities of two regions
  326. def CloseLoopSensitivities_TwoRegions(Sn_fp,Sg_fp,autocoupling_E,coupling_ItoE_abs,basal_forcing_E,coupling_EtoI,autocoupling_I_abs,basal_forcing_I,Tglu,Tgaba,G,Jnmda,k,M,control,target):
  327. Kappa = G*Jnmda*M
  328. Currents_OtherRegions_ToE = dot(k*Kappa,Sn_fp)
  329. k_inh = get_k_inh(k)
  330. Currents_OtherRegions_ToI = dot(k_inh*Kappa,Sn_fp)
  331. #--- Region 1 (open loop)
  332. i=control
  333. (SnC,SgC,coupling_ItoE_abs_C,TgabaC,TgluC) = get_ParametersForRegion_i(i,Sn_fp,Sg_fp,coupling_ItoE_abs,Tgaba,Tglu)
  334. xEC = basal_forcing_E[i] + Currents_OtherRegions_ToE[i]
  335. xIC = basal_forcing_I + Currents_OtherRegions_ToI[i]
  336. (Sensi_SnC_xEC,Sensi_SgC_xEC,Sensi_SnC_xIC,Sensi_SgC_xIC) = OpenLoopSensitivities_OneRegion(SnC,SgC,xEC,xIC,autocoupling_E,coupling_ItoE_abs_C,coupling_EtoI,autocoupling_I_abs,TgluC,TgabaC)
  337. #--- Region 2 (open loop)
  338. i=target
  339. (SnT,SgT,coupling_ItoE_abs_T,TgabaT,TgluT) = get_ParametersForRegion_i(i,Sn_fp,Sg_fp,coupling_ItoE_abs,Tgaba,Tglu)
  340. xET = basal_forcing_E[i] + Currents_OtherRegions_ToE[i]
  341. xIT = basal_forcing_I + Currents_OtherRegions_ToI[i]
  342. (Sensi_SnT_xET,Sensi_SgT_xET,Sensi_SnT_xIT,Sensi_SgT_xIT) = OpenLoopSensitivities_OneRegion(SnT,SgT,xET,xIT,autocoupling_E,coupling_ItoE_abs_T,coupling_EtoI,autocoupling_I_abs,TgluT,TgabaT)
  343. #--- Feedback gain
  344. (G12,G21,Gain) = FeedbackGain_TwoRegions(k,Kappa,Sensi_SnC_xEC,Sensi_SnC_xIC,Sensi_SnT_xET,Sensi_SnT_xIT)
  345. #--- Close loop system
  346. (dSnC_dBEC,dSnT_dBEC) = SensitivitiesTwoRegions_Sn_wrt_BEC(Sensi_SnC_xEC,Gain,G21)
  347. return(dSnC_dBEC, dSnT_dBEC, Gain, Sensi_SnC_xEC, Sensi_SgC_xEC, Sensi_SnC_xIC, Sensi_SgC_xIC, Sensi_SnT_xET, Sensi_SgT_xET, Sensi_SnT_xIT, Sensi_SgT_xIT)
  348. def get_ParametersForRegion_i(i,Sn,Sg,b,Tgaba,Tglu):
  349. Sni = Sn[i]
  350. Sgi = Sg[i]
  351. bi = b[i]
  352. Tgabai = Tgaba[i]
  353. Tglui = Tglu[i]
  354. return(Sni,Sgi,bi,Tgabai,Tglui)
  355. def FeedbackGain_TwoRegions(k,Kappa,Sensi_Sn1_xE1,Sensi_Sn1_xI1,Sensi_Sn2_xE2,Sensi_Sn2_xI2):
  356. G12 = Sensi_Sn1_xE1 * k[0,1]* Kappa[0,1] + Sensi_Sn1_xI1 * (1-k[0,1])* Kappa[0,1]
  357. G21 = Sensi_Sn2_xE2 * k[1,0] * Kappa[1,0] + Sensi_Sn2_xI2 * (1-k[1,0]) * Kappa[1,0]
  358. Gain = G12 * G21
  359. return(G12,G21,Gain)
  360. def SensitivitiesTwoRegions_Sn_wrt_BEC(Sensi_Sn1_xE1,Gain,G21):
  361. dSn1_dBE1 = (Sensi_Sn1_xE1)/(1-Gain)
  362. dSn2_dBE1 = G21 * dSn1_dBE1
  363. return(dSn1_dBE1,dSn2_dBE1)

functions.py, no license · at the source

Overview

Authors: Anaïs Vallet1, Stéphane Blanco2, Coline Chevallier2,3, Francis Eustache1, Jacques Gautrais2,3, Jean-Yves Grandpeix4, Jean-Louis Joly2, Shailendra Segobin1, Pierre Gagnepain1
  1. Normandie Univ, UNICAEN, PSL Research University, EPHE, INSERM, U1077, CHU de Caen, GIP Cyceron, Neuropsychologie et Imagerie de la Mémoire Humaine, Université de Normandie, Caen, France
  2. CNRS, INPT, UPS, LAPLACE, Université de Toulouse, Toulouse, France
  3. Centre de Recherches sur la Cognition Animale (CRCA), Centre de Biologie Intégrative (CBI), CNRS, UPS, Université de Toulouse, Toulouse, France
  4. LMD/IPSL, CNRS, École Polytechnique, ENS, Sorbonne Université, Paris, France
Journal: PLoS computational biology, volume 22, issue 3, article e1014035
Dates: received 30 April 2025; accepted 17 February 2026; published online 23 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014035 · PMID 41871132 · PMCID PMC13008111 · OpenAlex W7140113130
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: human (organism), computational (subfield)
MeSH: Brain*, Models, Neurological*, Nerve Net*, Computational Biology, Computer Simulation, Humans, Neural Inhibition, Neurons (* major topic)
Journal subjects: Biology and Life Sciences, Neuroscience, Brain Mapping, Functional Magnetic Resonance Imaging, Medicine and Health Sciences, Diagnostic Medicine, Diagnostic Radiology, Magnetic Resonance Imaging, Research and Analysis Methods, Imaging Techniques, Radiology and Imaging, Neuroimaging, Biophysics, Physical Sciences, Physics, Cell Biology, Cellular Types, Animal Cells, Neurons, Cellular Neuroscience, Computer and Information Sciences, Neural Networks, Psychology, Behavior, Social Sciences, Signal Transduction, Cell Signaling, Signal Inhibition, Engineering and Technology, Signal Processing, Peak Values, Signal Amplification
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Citations: not cited yet (Europe PMC); 36 references in the paper

Abstract

Biophysical models of human brain represent the latter as a graph of inter-connected neural regions. Building from the model by Naskar et al. (Network Neuroscience, 2021), our motivation was to understand how these brain regions can be connected at neural level to implement some inhibitory control, which calls for inhibitory connectivity rarely considered in such models. In this model, regions are made of inter-connected excitatory and inhibitory pools of neurons, but are long-range connected only via excitatory pools (mutual excitation). We thus extend this model by generalizing connectivity, and we analyse how connectivity affects the behaviour of this model. Focusing on the simplest paradigm made of a Control area and a Target area, we explore four typical kinds of connectivity: mutual excitation, Target inhibition by Control, Control inhibition by Target, and mutual inhibition. For this, we build an analytical sensitivity framework, nesting up sensitivities of isolated pools, of isolated regions, and of the full system. We show that inhibitory control can emerge only in Target inhibition by Control and mutual inhibition connectivities. We next offer an analysis of how the model sensitivities depends on connectivity structure, depending on a parameter controling the strength of the self-inhibition within Target region. Finally, we illustrate the effect of connectivity structure upon control effectivity in response to an external forcing in the Control area. Beyond the case explored here, our methodology to build analytical sensitivities by nesting up levels (pool, region, system) lays the groundwork for expressing nested sensitivities for more complex network configurations, either for this model or any other one.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above.

supp:PMC13008111/pcbi.1014035.s003.zip

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Languages: Python (11), Shell (1)
Size: 12 files, 12 scripts
Software Heritage: not checked
Found in: the supplementary material
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (11 files), Matplotlib (10 files)
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
12 files

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 12 scripts, each with its path and the digest of its content;
  • no match between paragraphs and code yet;
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

All relevant data are within the manuscript and its Supporting Information files.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 8 MeSH terms, 3 funders, 36 references.

Cite

This paper

Vallet, A., Blanco, S., Chevallier, C., Eustache, F., Gautrais, J., Grandpeix, J.-Y., Joly, J.-L., Segobin, S., & Gagnepain, P. (2026). Sensitivity analysis enlightens effects of connectivity in a Neural Mass Model under Control-Target mode. PLoS computational biology, 22(3), e1014035. https://doi.org/10.1371/journal.pcbi.1014035

BibTeX

@article{vallet2026sensitivity,
author = {Vallet, Anaïs and Blanco, Stéphane and Chevallier, Coline and Eustache, Francis and Gautrais, Jacques and Grandpeix, Jean-Yves and Joly, Jean-Louis and Segobin, Shailendra and Gagnepain, Pierre},
title = {{Sensitivity analysis enlightens effects of connectivity in a Neural Mass Model under Control-Target mode}},
journal = {PLoS computational biology},
year = {2026},
month = mar,
volume = {22},
number = {3},
pages = {e1014035},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014035},
url = {https://doi.org/10.1371/journal.pcbi.1014035},
pmid = {41871132},
pmcid = {PMC13008111}
}

RIS

TY - JOUR
AU - Vallet, Anaïs
AU - Blanco, Stéphane
AU - Chevallier, Coline
AU - Eustache, Francis
AU - Gautrais, Jacques
AU - Grandpeix, Jean-Yves
AU - Joly, Jean-Louis
AU - Segobin, Shailendra
AU - Gagnepain, Pierre
TI - Sensitivity analysis enlightens effects of connectivity in a Neural Mass Model under Control-Target mode
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/03/23
VL - 22
IS - 3
SP - e1014035
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014035
UR - https://doi.org/10.1371/journal.pcbi.1014035
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014035",
"type": "article-journal",
"title": "Sensitivity analysis enlightens effects of connectivity in a Neural Mass Model under Control-Target mode",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Vallet",
"given": "Anaïs"
},
{
"family": "Blanco",
"given": "Stéphane"
},
{
"family": "Chevallier",
"given": "Coline"
},
{
"family": "Eustache",
"given": "Francis"
},
{
"family": "Gautrais",
"given": "Jacques"
},
{
"family": "Grandpeix",
"given": "Jean-Yves"
},
{
"family": "Joly",
"given": "Jean-Louis"
},
{
"family": "Segobin",
"given": "Shailendra"
},
{
"family": "Gagnepain",
"given": "Pierre"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "3",
"page": "e1014035",
"DOI": "10.1371/journal.pcbi.1014035",
"PMID": "41871132",
"PMCID": "PMC13008111",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014035",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
23
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.isci.2022.105516 [code]
Attentional capture mediates the emergence and suppression of intrusive memories.
Journal: iScience
In common: 5 references, 2 authors
[2] doi:10.1126/sciadv.adq8336 [code]
Plasticity of human resilience mechanisms.
Journal: Science advances
In common: 4 references, 2 authors
[3] doi:10.1038/s41467-022-30855-x [code]
Altered predictive control during memory suppression in PTSD.
Journal: Nature communications
In common: computational, 3 references, 2 authors
[4] doi:10.7554/elife.106194 [code]
Virtual Brain Inference (VBI), a flexible and integrative toolkit for efficient probabilistic inference on whole-brain models.
Journal: eLife
In common: Matplotlib, NumPy, computational, 10 references
[5] doi:10.1002/advs.202406440 [code]
TVB C++: A Fast and Flexible Back-End for The Virtual Brain.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: Matplotlib, NumPy, 8 references
[6] doi:10.1038/s41598-024-83542-w [code]
Symmetry breaking organizes the brain's resting state manifold.
Journal: Scientific reports
In common: computational, 8 references
[7] doi:10.1371/journal.pcbi.1012595 [code]
Homeodynamic feedback inhibition control in whole-brain simulations.
Journal: PLoS computational biology
In common: Matplotlib, NumPy, computational, 7 references
[8] doi:10.1016/j.neuroimage.2023.119926 [code]
Reduced emergent character of neural dynamics in patients with a disrupted connectome.
Journal: NeuroImage
In common: NumPy, 8 references
[9] doi:10.1371/journal.pcbi.1011434 [code]
A multi-layer mean-field model of the cerebellum embedding microstructure and population-specific dynamics.
Journal: PLoS computational biology
In common: Matplotlib, NumPy, computational, 6 references
[10] doi:10.1162/netn_a_00410 [code]
Neural mass modeling for the masses: Democratizing access to whole-brain biophysical modeling with FastDMF.
Journal: Network neuroscience (Cambridge, Mass.)
In common: NumPy, computational, 7 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1371/journal.pcbi.1014035.svg)](https://openscicode.org/badge/doi/10.1371/journal.pcbi.1014035)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1371/journal.pcbi.1014035.svg
   :target: https://openscicode.org/badge/doi/10.1371/journal.pcbi.1014035
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1371/journal.pcbi.1014035"><img src="https://openscicode.org/badge/doi/10.1371/journal.pcbi.1014035.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.