使用PyRates进行基于Jansen-Rit模型的神经动力学建模与仿真
PyRates是一个基于神经元放电率的神经动力学仿真库，使用该库可以实现基于Jansen-Rit模型的神经动力学模型建模。
PyRates的安装
PyRates可以使用`pip`命令进行安装：
pip install pyrates在Python代码文件中，引入pyrates包即可使用PyRates：
import pyrates
基本结构
Jansen-Rit模型环路（Circuit）中的每个环路节点包含3个子节点：
- 兴奋性神经元（excitatory interneurons, EIN）
- 抑制性神经元（inhibitory interneurons, IIN）
- 锥体细胞（Pyramidal cells, PC）
在PyRates中，每个子节点被描述为一系列“算子（Operator）”的集合，包括：
- 放电率-电位算子（Rate-to-potential operator, RPO）：将放电率转化为可测量的组织电位
- 电位-放电率算子（Potential-to-rate operator, PRO）：将组织电位转换为放电率
通常，EIN和IIN只包含RPO和PRO各一个，而PC则包含2个RPO（分别接受来自EIN和IIN的输出）和1个PRO。
EIN、IIN和PC构成一个环路节点，环路节点之间通过“边（Edge）”连接。
因此，在一个Jansen-Rit模型中，存在以下信息流（假设输入到PC的外部输入为Vi）：
- EIN：(前级输出+PC输出) → EIN兴奋性RPO → 电压 → PRO → PC兴奋性RPO → V1
- IIN：(前级输出+PC输出) → IIN兴奋性RPO → 电压 → PRO → PC抑制性RPO → V2
- PC：(V1+V2+Vi) → 电压 → PRO → PC输出（后级输入）
定义算子
PyRates将算子描述为pyrates.frontend.OperatorTemplate对象的实例，该对象通过equations属性定义其内部的数学方程组，并使用variables属性定义方程组中变量或常量的含义（如输入量"input"、输出量"output"、临时变量"variable"或具有固定数值的常量等，其中输出量有且仅能有1个）。
定义变量时，不能使用某些名称或后缀，因为它们为PyRates内部保留或预定义的符号，该非法符号列表由pyrates.frontend.OperatorTemplate.check_vname()函数定义：
disallowed_names = ['y', 'dy', 'source_idx', 'target_idx', 'pi', 'I', 'E'] disallowed_name_parts = ['_buffer', '_delays', '_maxdelay', '_idx', '_hist']
下面的方法可以用来定义兴奋性RPO（RPO_e）、抑制性RPO（RPO_i）和PRO的模板对象，从而避免反复定义算子方程组和变量：
from pyrates.frontend import OperatorTemplate, NodeTemplate, EdgeTemplate, CircuitTemplate # String constants # Operators sRPEDescTemplate = "Excitatory rate-to-potential operator" sRPIDescTemplate = "Inhibitory rate-to-potential operator" sPRODescTemplate = "Sigmoidal potential-to-rate operator" # Sub nodes sEINDescTemplate = "Excitatory interneurons (EIN)" sIINDescTemplate = "Inhibitory interneurons (IIN)" sPMCDescTemplate = "Pyramidal cells (PC)" # Jansen-Rit circuits sJRCDescTemplate = "Jansen-Rit circuit model" # Operator templates # Use copy.deepcopy() to create replica of these operator templates for your own nodes rpeTemplate = OperatorTemplate(name="RPO_e", path=None, equations=[ "d/dt * V = I_c", "d/dt * I_c = H/tau * m_in - 2 * I_c/tau - V/(tau^2)" ], variables={ "V" : "output", "I_c" : "variable", "m_in" : "input", "tau" : 0.01, "H" : 0.00325, }, description=f"{sRPEDescTemplate} template") rpiTemplate = OperatorTemplate(name="RPO_i", path=None, equations=[ "d/dt * V = I_c", "d/dt * I_c = H/tau * m_in - 2 * I_c/tau - V/(tau^2)" ], variables={ "V" : "output", "I_c" : "variable", "m_in" : "input", "tau" : 0.02, "H" : -0.022, }, description=f"{sRPIDescTemplate} template") proTemplate = OperatorTemplate(name="PRO", path=None, equations=["m_out = 2.*m_max / (1 + exp(r * (V_thr - V)))"], variables={ "m_out" : "output", "V" : "input", "V_thr" : 6e-3, "m_max" : 2.5, "r" : 560.0 }, description=f"{sPRODescTemplate} template")
使用copy.deepcopy()建立算子实例的深拷贝对象，PyRates的Template类提供了用于修改对象参数的update_template()接口：
import copy sCurrentNodeName = "Demo Jansen-Rit Circuit (DemoJRC)" # Operators rpeDemoJRC = copy.deepcopy(rpeTemplate).update_template(name="RPO_e", path=None, description=f"{sRPEDescTemplate} of {sCurrentNodeName}") rpiDemoJRC = copy.deepcopy(rpiTemplate).update_template(name="RPO_i", path=None, description=f"{sRPIDescTemplate} of {sCurrentNodeName}") proDemoJRC = copy.deepcopy(proTemplate).update_template(name="PRO", path=None, description=f"{sPRODescTemplate} of {sCurrentNodeName}")
定义子节点
PyRates将子节点描述为pyrates.frontend.NodeTemplate对象的实例，该对象通过operators属性定义其使用的算子列表）。
# Sub nodes einDemoJRC = NodeTemplate(name="EIN", path=None, operators=[rpeDemoJRC, proDemoJRC], description=f"{sEINDescTemplate} of {sCurrentNodeName}") iinDemoJRC = NodeTemplate(name="IIN", path=None, operators=[rpeDemoJRC, proDemoJRC], description=f"{sIINDescTemplate} of {sCurrentNodeName}") pmcDemoJRC = NodeTemplate(name="PC", path=None, operators=[rpeDemoJRC, rpiDemoJRC, proDemoJRC], description=f"{sPMCDescTemplate} of {sCurrentNodeName}")
定义环路节点
PyRates将环路节点描述为pyrates.frontend.CircuitTemplate对象的实例，该对象使用nodes属性定义节点中的子节点，使用circuits属性定义子环路，并使用edges属性定义环路中各边的起止点和权重等参数。
其中，nodes和circuits属性使用键值对的形式，定义子节点名称和对象实例之间的对应关系。edges属性使用pyrates.frontend.CircuitTemplate对象中子节点及算子的name和variables属性，确定该边从何输出变量获取数据，并传送到何输入变量。
# Sub node interconnections in Jansen-Rit Model jrcDemoJRC = CircuitTemplate(name="JRC_DemoJRC", path=None, nodes={ "EIN" : einDemoJRC, "IIN" : iinDemoJRC, "PC" : pmcDemoJRC }, edges=[ # EIN -> PC ("EIN/PRO/m_out", "PC/RPO_e/m_in", None, {"weight" : 108}), # IIN -> PC ("IIN/PRO/m_out", "PC/RPO_i/m_in", None, {"weight" : 33.75}), # PC -> EIN ("PC/PRO/m_out", "EIN/RPO_e/m_in", None, {"weight" : 135}), # PC -> IIN ("PC/PRO/m_out", "IIN/RPO_e/m_in", None, {"weight" : 33.75}), ], description=f"{sJRCDescTemplate} of {sCurrentNodeName}")
开始仿真
使用pyrates.frontend.CircuitTemplate对象的run()函数开始仿真，该函数的常用参数包括：
- simulation_time：仿真时长，单位为秒。
- step_size：仿真步长，单位为秒。
- sapling_step_size：采样步长，单位为秒。
- outputs：仿真输出键值对，键名为仿真输出（虚拟采样通道）的索引名称，值为pyrates.frontend.CircuitTemplate对象中子节点及算子的name和variables属性确定的输出名称，表示记录来自哪个变量的数据。注意：只有输出变量（类型为“output”的变量）可被测量。
- inputs：输入键值对，键名为pyrates.frontend.CircuitTemplate对象中子节点及算子的name和variables属性确定的输入名称，表示向哪个变量输入数据。值为预先生成的输入信号（一维数组），输入信号的采样间隔应与仿真步长一致（即输入信号数组的长度应为仿真时长除以仿真步长）。
dSimulationDuration = 5.0 dSimulationStep = 1e-4 dSamplingStep = 1.0 / 500.0 dSimulationPaddingAtStart = 0.0 # Calculate real simulation duration # PyRates outputs simulation_time*(1.0/sampling_step_size) points after simulation # Thus the real timespan (MNE algorithm) is simulation_time-dSamplingStep seconds # To fix this, you can add 1 sampling_step_size to simulation_time. # dSimulationPaddingAtStart is adopted to contain unstable leading segment dSimulationDurationReal = dSimulationPaddingAtStart + dSimulationDuration + dSamplingStep # Generate input signal (square wave) # Square wave is generated using scipy.signal.square() # Note that scipy.signal.square() generates a square wave with a period of 2*pi and a range of [-1, 1] # See https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.square.html#scipy.signal.square # In this demo, we generated a signal with 130 Hz frequency and 10% cycle duty, in the range of [0, 5] dInputSignalAmplitude=2.5 dInputSignalFreq=130 dInputSignalPhase=0.0 dInputSignalDuty=0.1 dInputSignalOffset=2.5 arrInputSignalX = np.linspace(0, dSimulationDurationReal, round(dSimulationDurationReal/dSimulationStep), endpoint=False) arrInputSignalY = scipy.signal.square((2 * np.pi * dInputSignalFreq * arrInputSignalX + dInputSignalPhase - dSimulationPaddingAtStart), duty=dInputSignalDuty) arrInputSignalY = dInputSignalAmplitude * arrInputSignalY + dInputSignalOffset # Run simulation # Results is an instance of pandas.DataFrame, whose columns are output channels, indexes are time points # Use the "inputs" parameter to specify input signal # Use the "outputs" parameter to specify output channels, you can access them with [ChannelName] index # Only "output" variables defined in OperatorTemplate can be measured # See https://pyrates.readthedocs.io/en/latest/frontend.template.html#pyrates.frontend.template.circuit.CircuitTemplate.run dctSimResults = jrcDemoJRC.run( simulation_time=dSimulationDurationReal, step_size=dSimulationStep, sampling_step_size=dSamplingStep, inputs={ # Stimulating PC, input to operator PRO, variable V "PC/PRO/V" : arrInputSignalY } outputs={ # Record from subnode PC, operator RPO_e (EIN output), variable V "PC_e" : "PC/RPO_e/V" # Record from subnode PC, operator RPO_i (IIN output), variable V "PC_i" : "PC/RPO_i/V" }, clear=True )
需要注意的是，PyRates生成的仿真结果中，共包含nPoint = simulation_time * (1.0 / sampling_step_size)个测量点。因此，部分算法（例如MNE工具包采用的算法：dTime = (nPoint - 1) * sampling_step_size）会导致测得的时长略小于simulation_time的设置值。此时，可以通过将simulation_time设为期望的仿真时长与1个sampling_step_size的和进行处理，即设置simulation_time=dSimulationDuration+dSamplingStep。
为了获得较为稳定的仿真结果，可以考虑在simulation_time中设置一个5 s左右的引导段，并在处理时剔除这个引导段：
# Remove simulation padding segment iSimulationValidDataStartIndex = dSimulationPaddingAtStart / dSamplingStep iSimulationValidDataStartIndex = round(dSimulationPaddingAtStart / dSamplingStep) arrSimulationDataIndexesToDrop = list(range(0, iSimulationValidDataStartIndex)) arrSimulationDataTimePointsToDrop = [i*dSamplingStep for i in arrSimulationDataIndexesToDrop] dctSimResults = dctSimResults.drop(labels=arrSimulationDataTimePointsToDrop, axis=0) dctSimResults = dctSimResults.reset_index(drop=True) arrSimualtionDataNewIndex = [i*dSamplingStep for i in dctSimResults.index] dctSimResults.index = arrSimualtionDataNewIndex
由于来自同一环路的EIN和IIN的电位被分开测量，因此可以使用加法对其进行合并：
# Merging results from EIN and IIN dctSimResultsMerged = pd.DataFrame() for CurrentKey in dctSimResults.keys() : # Check if current key is from a "XXX_e" and "XXX_i" pair if str(CurrentKey).endswith("_e") : # If curren key is an EIN output sCurrentKeyBase = str(CurrentKey).removesuffix("_e") # Dertermin key of paired IIN output sCurrentKeyPaired = sCurrentKeyBase + "_i" if sCurrentKeyPaired in dctSimResults.keys() : # Merging key pairs, a simple "plus" operator is used # See https://pyrates.readthedocs.io/en/latest/auto_implementations/python_definitions.html dctSimResultsMerged[sCurrentKeyBase] = dctSimResults[sCurrentKeyBase + "_e"] + dctSimResults[sCurrentKeyBase + "_i"] else : # If we can't find a paired IIN output, log this entry directly dctSimResultsMerged[CurrentKey] = dctSimResults[CurrentKey] #End If elif str(CurrentKey).endswith("_i") : # If curren key is an IIN output sCurrentKeyBase = str(CurrentKey).removesuffix("_i") # Dertermin key of paired EIN output sCurrentKeyPaired = sCurrentKeyBase + "_e" if sCurrentKeyPaired in dctSimResults.keys() : # Since EIN output is handled in previous If branch, simply continuing the loop continue else : # If we can't find a paired EIN output, log this entry directly dctSimResultsMerged[CurrentKey] = dctSimResults[CurrentKey] #End If else : # For entries with other suffixes, log this entry directly dctSimResultsMerged[CurrentKey] = dctSimResults[CurrentKey] #End If #Next
仿真结果被存储为pandas.Dataframe对象的实例，其数据列（Columns）为输出键名，数据索引（Indexes）为数据对应的时间点。您可以通过run()函数的outputs中指定的键名访问仿真过程中获得的采样数据。
import matplotlib.pyplot as plt import pandas as pd # Plotting data for CurrentChannel in dctSimResultsMerged.keys() : plt.plot(dctSimResultsMerged[CurrentChannel]) plt.title(CurrentChannel) #Next # Exporting data dctSimResultsMerged.to_csv("SimResults.csv", index=False)
杂项
缓存文件问题
PyRates默认会在当前工作目录（通过os.getcwd()获取）下建立缓存文件（如pyrates_run.py），若需要并行执行多个PyRates仿真，可能需要动态调整工作目录，例如：
# Change working directory to avoid pyrates_run.py confliction import os import shutil import uuid sSimulatorCacheDir = "temp/sim/" sPreviousWorkDir = os.getcwd() sCurrentSimulatorCacheDir = sSimulatorCacheDir + str(uuid.uuid4()) + "/" os.makedirs(sCurrentSimulatorCacheDir, exist_ok=True) os.chdir(sCurrentSimulatorCacheDir) # Running simulation jrcModel.run(...) # Revert working directory os.chdir(sPreviousWorkDir) shutil.rmtree(sCurrentSimulatorCacheDir, ignore_errors=True)
参考资料
https://pyrates.readthedocs.io/en/latest/auto_implementations/python_definitions.html
https://pyrates.readthedocs.io/en/latest/auto_introductions/jansenrit.html
页面版本: 23, 最后编辑于: 06 Dec 2025 02:42