Pigreads：一款集成Python语言的、基于GPU的反应-扩散求解器，支持OpenCL技术，广泛应用于心脏电生理学及其他相关领域

《Computer Physics Communications》：Pigreads: The Python-integrated GPU-enabled reaction-diffusion solver using OpenCL for cardiac electrophysiology and other applications

　　德斯蒙德·卡布斯（Desmond Kabus）|汉斯·迪尔克克斯（Hans Dierckx）|蒂姆·德科斯特（Tim De Coster）
莱顿大学医学中心（LUMC）实验心脏病学实验室，阿尔比努斯德雷夫2号（Albinusdreef 2），莱顿，2333 ZA，荷兰

**摘要**
Pigreads是一个高效的Python模块，用于在图形卡（GPU）上数值求解反应-扩散系统，并具有CPU回退功能。它提供了一个简单直接的、兼容NumPy的API。用户可以使用内置模型（包括电生理学示例），或提供自定义的反应项。支持的功能包括0D-3D均匀笛卡尔网格、无通量边界条件、各向异性扩散、空间变化的扩散和反应以及局部源项。该项目是开源的，已经过测试，并附带示例和教程。

**程序概述**
- **程序名称：** Pigreads
- **程序文件链接：** https://doi.org/10.17632/5k8jjx74hc.1
- **开发者仓库链接：** https://gitlab.com/pigreads/pigreads
- **许可条款：** MIT
- **编程语言：** Python, OpenCL
- **问题性质：** 反应-扩散系统用于模拟物理、化学和生物学中的现象；其数值求解在计算上可能非常耗时，且研究代码通常都是临时的、文档记录不足的。
- **解决方法：** 该实现采用有限差分空间离散化和显式前向欧拉时间步进法。性能关键的核心代码在OpenCL中运行，并通过PyOpenCL从Python调用；数据用NumPy数组表示。

**附加说明（包括限制和不寻常的功能）：**
Pigreads主要针对心脏电生理学应用的设计，强调简单性和可重复性，而非高级的数值复杂性。我们遵循科学软件开发的最佳实践，包括版本控制、测试和持续集成。

**1. 引言**
反应-扩散系统广泛应用于物理、化学和生物学领域，用于研究模式形成和兴奋波等现象（例如在电生理学中[1]、[2]、[3]、[4]、[5]、[6]）。其数值求解在计算上可能非常耗时，因此高效的实现至关重要，从临时的研究代码到大型通用软件包都有（[7]）。现代应用（如反问题、参数估计、不确定性量化和机器学习）需要一个占用空间小、安装要求低、跨平台可移植、运行速度快且易于集成到其他软件或框架中的代码。Pigreads是一个Python模块[8]，提供了一个简洁且兼容NumPy的API[9]，用于定义和运行反应-扩散模拟，并在加速器上高效执行计算：用OpenCL[10]实现的计算核心既可以在GPU上运行，也可以通过PyOpenCL在CPU上运行作为回退选项。Pigreads模块包含几个与电生理学相关的预定义模型；用户还可以提供自定义的反应项。
**代表性示例：** 图1展示了在一维电缆中的传导模拟、培养孔中的二维心脏组织片以及三维心脏器官（如心室或心房）的数值电生理实验。也支持无扩散的情况（0D情形），以进行单细胞模拟。

**下载：**
- 下载高分辨率图片（2MB）
- 下载全尺寸图片

**图1.** Pigreads可用于求解最多三个维度的反应-扩散系统。
A. Courtemanche、Ramirez和Nattel模型的0D情形下的心脏动作电位。
B. 电缆中的三个行波脉冲，随后是同一个模型中的第四个阻滞脉冲。
C. Marcotte和Grigoriev模型中的圆形二维域中的螺旋波，通过改变恢复变量的高低值来刺激。
D. ten Tusscher和Panfilov模型中的二维组织片中的伪电图，顶部和底部具有周期性边界条件，中间壁、内膜和心外膜细胞类型不同，并在标记为S0–2的位置和时间点刺激了三个脉冲。
E. 三维心脏电生理学模拟，分别在人类双心室[15]、[16]和双心房[18]几何结构中，肌肉纤维中的扩散速度更快。从心尖或窦房结分别发送一个脉冲。我们使用ten Tusscher和Panfilov[14]的模型来模拟心室，以及Courtemanche、Ramirez和Nattel[12]的模型来模拟心房。

该模块可在[https://gitlab.com/pigreads/pigreads]获取，已通过持续集成进行了充分测试，并在[Read the Docs (https://pigreads.readthedocs.io)]中提供了详细文档和示例教程。可以通过Python包索引PyPI（https://pypi.org/project/pigreads）安装：
```bash
pip install pigreads
```
**还建议安装可选的命令行接口（CLI）依赖项，以便绘制图表和视频：**
```bash
pip install pigreads[all]
```
CLI可用于运行和可视化根据Pydantic（https://github.com/pydantic/pydantic）规范定义在YAML文件[20]中的Pigreads模拟。

**2. 数学问题定义**
我们通常将反应-扩散系统定义为：
对于时间t ∈ [0, T]和空间x ∈ ?3中的域?：
$$
\partial_t u'(t, x) = P \cdot \nabla \cdot D(x) \nabla u'(t, x) + r(u', x) + s(t, x)
$$
初始条件为u(0, x)，并且在x ∈ ??上的扩散变量满足无通量边界条件0 = n \cdot D \nabla u"。
在每个时间和空间点，状态向量u'(t, x) ∈ ?^n包含Nv个元素（即状态或变量u0, u1, ..., uNv?1）。对于扩散项，定义扩散矩阵D(x) ∈ ?3×3和选择矩阵P ∈ ?^n×Nv，用于选择哪些变量进行扩散以及扩散强度。选择矩阵通常采取稀疏形式P = diag(P0, ..., PNv?1)。反应项r(u', x) ∈ ?^n描述系统的局部动态，可能在空间上变化（例如在参数或模型方程本身中）。我们将具有固定参数的特定局部选择r'imodel(u')称为一个模型；r表示所有模型。源项s(t, x) ∈ ?^n可用于向系统添加外部影响，例如在特定时间和位置刺激系统。
在没有源项（s = 0）的情况下，为均匀各向同性扩散，D(x) = const，且只有两个变量u = u0和v = u1（只有u扩散），则系统简化为：
$$
\partial_t u(t, x) = D \nabla^2 u(t, x) + ru(u, v)
$$

**3. 技术概述**
Pigreads围绕三个对象组织模拟组件：
- `Models`类实例，用于定义局部反应项；
- `inhom`数组，用于定义域和使用的模型；
- `states`数组，用于存储均匀笛卡尔网格上的所有状态变量；
**参见图2的流程图。**

**4. 动手示例**
在本节中，我们定义并运行一个简单的二维模拟，模拟在具有无通量边界条件的圆形域中的螺旋波。通过这个例子，可以了解通常如何定义Pigreads模拟的各个步骤：
首先，定义要使用的几何坐标。在这个例子中，我们使用一个包含200个点（x和y方向各100个点）的二维平面：
```python
import pigreads as pig
import numpy as np

R = 10
z, y, x = np.mgrid[0:1, -R:R: 200]
Ny, Nx = x.shape
z, dy, dx = pig.deltas(z, y, x)
```
Pigreads针对三维空间进行了优化。对于低维模拟，将其他维度的点数设置为1（如z维度所示）。注意，`np.mgrid`用于定义一个密集的多维网格，其中坐标范围是从0到R（不包括1），并且在y和x方向上分成200个等份。
Pigreads中的空间默认是周期性的，即索引ix = 0的点与ix = 1和ix = Nx - 1的点是邻居。
网格间距dz、dy、dx必须选择得足够小以便解析系统的动态，同时又足够大以保持计算成本可控。例如，如果空间分辨率不足，沿网格轴的传导速度会比对角线慢。可以利用这一效应来确定尽可能大的网格间距dx，使得在某个距离内最快和最慢的刺激之间的差异小于给定阈值（例如5%）。
`inhom`整型字段用于指定哪些点在域?之外（inhom = 0）或之内（inhom > 0）。Pigreads在??上实现无通量边界条件。在我们的例子中，我们将域?定义为半径为R的圆盘：
```python
inhom = np.ones(x.shape, dtype=int)
r = np.linalg.norm((x, y, z), axis=0)
inhom[r >= R] = 0
```
大于0的inhom值可用于选择一个或多个模型（即反应项r）。例如，inhom = 1时使用models[0]模型；inhom = 2时使用models[1]模型等。可以使用`Models`类的实例选择一个或多个模型：
```python
models = pig.models()
models.add('marcotte2017dynamical', beta=1.389)
Nv = models.Nv
```
模型的键（一个标识字符串）用于选择模型，关键字参数用于将模型参数设置为默认值之外的值；可以通过多次调用`add`函数来添加更多模型。当使用多个模型时，将使用最多的变量数量Nv。
Pigreads中预定义了多种模型，定义新模型非常简单（见第5节）。还可以通过以下代码程序获取可用模型列表：
```python
for key in pig.models.available.keys():
print(key)
```
**例如：**
```python
# aliev1996simple
# barkley1991model
# beeler1977reconstruction
# ...
```
对于某些模型，有多组参数可供使用，可以通过模型的元数据访问：
```python
key = 'tentusscher2006alternans'
model_def = pig.models.available[key]
endo = model_def.meta['parameter sets']['endo']
for parameter, value in endo.items():
print(f'{parameter}={value}')
# 示例参数值
# g_Ks = 0.392
# g_to = 0.073
# s_offset = 28.0
# s_variant = 1.0
models_ = pig.models()
models_.add(key, **endo)
```
接下来，需要分配内存来存储模型的状态变量。在Pigreads中，通常使用一个形状为(Nfr, Nz, Ny, Nx, Nv)的5D数组`states`来存储这些变量。
一种初始化`states`数组的方法是使用`models.resting_states`函数：它根据`inhom`的值为每个模型创建正确形状的数组，并填充第一个帧的初始值：
```python
Nfr = 100
states = models.resting_states(inhom, Nframes=Nfr)
assert states.shape == (Nfr, Nz, Ny, Nx, Nv)
```
需要注意的是，这会分配大量内存（因为`states`是一个形状为(Nfr, Nz, Ny, Nx, Nv)的5D数组）。对于大规模仿真，可能需要只存储少量的帧并覆盖较旧的帧，或者将帧保存到磁盘而不是保留在内存中，例如通过内存映射来实现：
```python
states = np.lib.format.open_memmap('states.npy', mode='w+', dtype=np.float32, shape=(Nfr, Nz, Ny, Nx, Nv)
states[:1] = models.resting_states(inhom, Nframes=1)
```
如果不需要保留完整的時間历史记录，也可以创建一个只包含少数帧的 `states` 数组，并在无限循环中不断覆盖它们。我们在代码仓库中的一些交互式示例中使用了这种方案。然后可以在第 `ifr` 帧中设置初始条件：
```python
ifr = 0:
states[0, x < -8, 0] = 1
states[0, y < 0, 1] = 2
```
或者用数学符号等价表达为：
```math
(4) u0(t=0, x) = {1 if x < -8 else u0
(5) u1(t=0, x) = {2 if y < 0 else u1
```
扩散项 `P?·D?u` 的计算是通过相邻点的加权和来实现的。权重可以使用 `weights` 函数计算，该函数还需要输入扩散系数 `D`：
```python
diffusivity = pig.diffusivity_matrix(Df=0.03)
weights = models_weights(dz, dy, dx, inhom, diffusivity)
```
计算权重是一个昂贵的操作。一旦为给定的几何形状计算出权重，就可以存储并重复使用它们。这些权重与 `inhom` 一起，完全编码了模型运行所需的全部几何信息。最后，可以使用 `run` 函数开始仿真，将仿真从一个帧推进到下一个帧，这在遍历帧数 `Nfr` 的循环中完成。`run` 函数执行 `Nt` 个前向欧拉步长，并在这些步骤之后返回最终状态：
```python
Nt = 200
dt = 0.025
for ifr in range(Nfr - 1):
states[ifr + 1] = models.run(inhom, weights, states[ifr], Nt=Nt, dt=dt)
```
注意，为了数值稳定性，时间步长 `dt` 需要选择得足够小，至少满足 Courant-Friedrichs-Lewy (CFL) 条件 [23]，非正式地可以总结为：“波在每个时间步长内传播的距离必须小于一个网格长度。” 更小的时间步长会带来更高的精度，但计算成本也会增加。我们通常选择尽可能大的 `dt`，同时仍然能够得到一个不会因 `dt` 减小而显著变化的解。`run` 函数的额外参数可用于在特定时间和位置添加刺激电流。有关更多详细信息，请参阅文档。

现在仿真已经完成，可以分析并可视化包含结果的 5D 数组 `states`，例如使用 `Matplotlib` [21]：
```python
import matplotlib.pyplot as plt
plt.imshow(states[-1, 0, :, :, 0])
plt.show()
```
或者使用 `FFmpeg` (https://ffmpeg.org) 将其作为电影显示：
```python
from pigreads.plot import moviemovie('example.mp4', states[:, 0, :, :, 0])
```
更多深入的示例可以在项目的综合 API 文档 (https://pigreads.readthedocs.io) 和 Git 仓库 (https://gitlab.com/pigreads/pigreads) 中找到。

5. 反应项
所谓的模型定义了反应-扩散方程 (Eq. (1)) 中的反应项 `r`。虽然 Pigreads 提供了多种预定义的模型，但也可以很容易地定义自定义模型。可以通过添加几行 OpenCL 代码将其添加到可用模型字典中来定义模型，例如 FitzHugh [24] 以及 Nagumo、Arimoto 和 Yoshizawa [25] 的模型：
```python
import pigreads as pig
from pigreads.schema.model import ModelDefinition
pig.models.available['fitzhugh1961impulses'] = \
ModelDefinition(
name='FitzHugh 1961 & Nagumo 1962',
description='A 2D simplification of the Hodgkin-Huxley model.',
doi=[['https://doi.org/10.1016/S0006-3495(61)86902-6', 'https://doi.org/10.1109/JRPROC.1962.288235'],
variables={'u': 1.2, 'v': -0.625},
diffusivity={'u': 1.0},
parameters={'a': 0.7, 'b': 0.8, 'c': 3.0, 'z': 0.0},
code=''''*_new_u = u + dt * (v + u - u*u*u/3 + z + _diffuse_u);
*_new_v = v + dt * (-(u - a + b*v)/c);
''''
)
```
要从 Myokit [26] 或 CellML [27] 文件导入模型，仓库中提供了相应的转换器。

虽然 Pigreads 的主要应用是心脏电生理学，但也有各种更通用的预定义模型：简单模型定义了没有反应项的单变量扩散；Gray 和 Scott [28] 的模型用于模式形成研究；Barkley [29] 以及 Hodgkin 和 Huxley [30] 的模型是一些最简单和最早的电生理学模型。这些模型的简单示例仿真如图 3 所示。

下载：下载高分辨率图片 (491KB)
下载：下载全尺寸图片

图 3. Pigreads 包含了多种预定义的模型。每个面板都包含了一个不同模型的 2D 仿真以及在域中标记的参考点处的时间迹线，这些参数是手动选择的。另见表 1. A. 简单模型定义了单变量的扩散。B. Gray 和 Scott [28] 的模型描述了根据所选参数形成的各种模式。C. Barkley [29] 的模型是一个最简单的模型，它能够产生螺旋波。D. Hodgkin 和 Huxley [30] 的模型描述了乌贼巨型神经纤维中的电传导，这是最早的电生理学数学模型之一。这里展示的是单个细胞的仿真，而不是时间迹线。

图 4 给出了心脏电生理学预定义模型的概览。这些仿真使用了自动确定的分辨率、时间和刺激，这些参数取决于模型的动态特性：时间步长 `dt`、网格间距 `dx`，以及通过瞬间将第一个变量增加给定值来施加的刺激幅度 `Δu`。在本节的其余部分，我们将概述这种自动寻优程序，以找到给定心脏电生理学模型的这些参数。用于此程序的 Python 脚本也在单独的第二个仓库 (https://gitlab.com/pigreads/figures) 中发布。

首先必须确定时间步长 `dt` 和刺激幅度 `Δu`。我们考虑一个可能值的范围 (dt, Δu) ∈ S ? [10^-8, 0) × [10^-5, 10^5]，并在其中运行单细胞仿真 (D=0)。然后我们将搜索范围 S 缩小到有效范围内仅包含一个动作电位的区域 S′ ? S。我们通过检查以下标准来确定：
- 所有变量在所有时间点上是否都是有限的？
- 第一个变量 `u` 是增加还是减少？
- `u` 的振荡次数是否少于五次？
- `u` 是否会偏离静息值？
- `u` 是否会在被刺激后超过激活值？
- `u` 的变化是否平滑，即没有大的跳跃？
- `u` 是否覆盖了多次仿真中观察到的大部分范围？

接着我们选择有效仿真 S′ 中第 80 百分位的 `dt`，然后在选定的 `dt` 下，选择有效仿真中第 30 百分位的 `Δu`。

然后我们运行另一个单细胞仿真来确定动作电位持续时间 (APD)，即 `u` 高于 30% 阈值的时间。

通过使用 CFL 条件 [23] 获得非常细的网格间距 `dx` 的 1D 仿真，我们可以确定传导速度 (CV)。在这里，我们暂时将扩散系数 `D` 设定为 1，因为之后可以按需要重新缩放它。

通过对第三种刺激在 1D 仿真中测量 CV 和 APD 来细化这些值。

在一系列 D=1 的 2D 仿真中，我们确定最大的网格间距 `dx`，使得沿轴线和对角线的最快速度和最慢速度之间的差异小于 5%。

如果在 2D 仿真中需要数值稳定性，即满足 CFL 条件，则为 `dt` 降低值。

图 4. Pigreads 中预定义的心脏电生理学模型。对于每个模型，我们展示了一个单个细胞在三脉冲之间的不同间隔下的跨膜电压仿真，以及使用 S0-S1-S2 协议启动的螺旋波的 2D 仿真。另见表 1. A.–E. 物理模型。F.–J. 心房细胞模型。K.–O. 心室细胞模型。注意：螺旋波在所有模型中不一定稳定，可能会分解成螺旋混沌。

这种自动程序的结果在表 1 中以表格形式给出了 Pigreads 中所有预定义模型的概览。该表包括了模型的原始出版物引用、变量数量 `Nv`、离散化参数 `dt` 和 `dx`、刺激强度 `Δu`，以及在扩散系数 D=1 时的 APD 和 CV。注意，这些参数可以通过以下公式缩放到所需的速度：
```math
D ∝ CV^2 / dx, dt ∝ const
```
这是根据方程 (1) 得出的，因为保持速度和扩散系数在网格单位中相同，即 D·dt/dx^2 和 CV·dt/dx 都是常数。

表 1. Pigreads 中预定义模型的概览。

| 作者 | 年份 | 参考文献 | Nv | dt (0D) | dt (2D) | dx | 刺激 | APD | CV |
|-----------| ----------|--------------|------|------|------|-----------|-----------------|---------------|
| Trivial model | * | 1982 | 10.00 | 40 | idem | 0.20 | 1.00 | |
| | | 1991 | 20.10 | idem | 1.00 | 0.40 | |
| | | 1991 | 20.00 | 20.0 | idem | 1.00 | 1.20 |
| | | 1952 | s | 40.00 | 50 | 15.0 | 2.25 |
| | | 1996 | p | 20.06 | 0.02 | 3.0 | 0.43 |
| | | 1998 | p | 30.85 | 0.21 | 1.3 | 25.4 |
| | | 2003 | p | 20.23 | idem | 3.3 | 0.08 |
| | | 2008 | p | 20.03 | 0.92 | 148.7 |
| | | 2008 | p | 20.03 | 0.92 | 148.7 |
| | | 2008 | p | 20.06 | 0.41 | 25.5 |
| | | 2017 | p | 20.03 | 0.92 | 148.7 |
| | | 2008 | p | 20.03 | 0.92 | 255.8 |
| | | 2008 | p | 20.06 | 0.41 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 148.7 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
| | | 2017 | p | 20.03 | 0.92 | 255.8 |
```

这种自动程序的结果在表 1 中以表格形式给出了 Pigreads 中所有预定义模型的概览。该模块也可用于教育目的，例如在数学生物学、动力系统或科学计算等课程中。虽然Pigreads可以通过多种方式进行扩展——例如采用更先进的时间步进方法、自适应网格细化、支持更复杂的几何形状或其他类型的控制方程——但我们选择保持模块的简洁性以尽可能确保其稳定性。虽然我们欢迎收到关于小修改、 bug 修复以及额外模型定义和几何文件的请求，但除非确实有必要，否则不应添加任何新功能。如果需要，始终可以分叉项目来添加额外功能，例如按项目具体需求进行扩展。我们希望Pigreads能成为各个科学领域研究人员的有用工具，并推动科学进步。

**资助**
本项工作得到了荷兰科学研究组织（NWO Open Mind 授予 TDC 的项目编号 2025/TTW/02025375）的支持。DK 得到了鲁汶大学（KU Leuven）的 GPUL/20/012 项目资助，以及莱顿大学的人工智能与生命科学研究所（SAILS）项目的支持。

**CRediT 作者贡献声明**
Desmond Kabus：撰写 – 评议与编辑、撰写 – 原始草稿、可视化、验证、软件开发、资源准备、方法学设计、数据分析、形式化分析、数据整理、概念构建。
Hans Dierckx：撰写 – 评议与编辑、指导工作。
Tim De Coster：撰写 – 评议与编辑、指导工作、项目管理工作、资金筹集、概念构建。