心脏,作为重要器官之一,其功能正常与否直接影响人类的生命延续。电生理特性反映了心脏的健康和疾病状态。心脏电生理活动的异常,往往会导致心律失常,从而引至心脏泵血功能衰竭等严重健康问题。因此,深入理解和研究心脏的电生理过程,对于提高心脏病的诊断和治疗水平至关重要。
传统的心脏电生理研究多依赖于实验室内的动物模型和临床数据,但这类方法往往受限于伦理问题、实验条件和数据获取的复杂性。随着计算技术的发展,计算机仿真成为了一种新兴且强大的研究工具。通过建立数学模型和计算机程序,研究人员可构建数字孪生心脏,能够在虚拟环境中仿真并重现心脏器官的电生理活动(虚拟生理心脏),分析其动态特性,并进行不同生理与病理条件下的实验。
虚拟心脏电生理仿真对计算资源要求极高,即使是几毫秒的仿真,也需要累积求解数十亿次微分方程。使用复杂的虚拟心脏模型进行研究时,重现1秒钟的心脏电活动也可能需要数小时或更长。这给虚拟生理心脏的临床应用与药物研发带来重大挑战。
为解决这一问题,智源研究院开发了一套实时心脏电生理仿真系统。该系统不仅能够实时模拟心脏的3D电活动,还能通过多种参数的调节,深入探讨不同生理、病理因素对心脏功能的影响。
这一实时心脏仿真平台,一方面可在医学基础研究领域发挥作用,帮助临床医生和研究人员更直观地理解心脏的电生理过程,探究心律失常产生机制、预测猝死发生率等;另一方面,可用于构建虚拟药物安全性评估平台,对推动药物安全评估发展具有重要意义;更重要的是,可以在临床应用中提供手术方案预演与决策支持,比如射频消融方案规划,心脏起搏器最佳植入方案规划等。该技术的推进将为医学研究和临床治疗提供新的范式。
1
虚拟心脏仿真发展史
虚拟生理心脏的构建可利用生理组学的研究方法,综合分子生物学、生物化学、生理学、解剖学及临床医学的最新成果,数学化以及模式化地整合从基因、蛋白质、细胞、组织到器官的解剖(多物理尺度:空间尺度跨越10^9量级,跨时间尺度:时间尺度跨越10^15量级,如图1所示)、生理和生化信息,应用计算机强有力的计算和图形显示能力,通过赋予其心脏所具有的动力学特性、生化特性和各种生理病理特点,使之从形态、结构和功能等方面逼真地再现心脏的生理和病理活动过程。
图注1:构建虚拟生化生理人体的时间和空间尺度。时间尺度横跨由分子事件(µs)、细胞信号传导(ms)、细胞功能(s)到人体寿命 (decades) 的10^15跨度。空间尺度横跨由分子(nm)、细胞(µm)、器官(cm)到躯干 (m) 的10^9跨越。
虚拟生理心脏研究可追溯于上世纪五十年代。1952年诺贝尔奖得主Hodgkin和Huxley建立了世界上第一个细胞计算模型 — 乌贼神经元细胞模型[1],开创了用计算模型研究生物问题的先河。1960年Denis Noble[2]在Nature杂志上发表了第一个心肌细胞计算模型 — 浦肯野心肌细胞模型,开创虚拟生理心脏模型的先例。此后几十年的研究中,不断有研究人员研发针对不同物种、心脏不同组织、复杂精密的心肌细胞电生理模型[3]。1991年,Peter Hunter等人[4]基于犬实验数据构建了第一个心脏解剖结构模型,融合多物理尺度与电生理的虚拟心脏模型研究进入新阶段。此后,多尺度、多物理模态的心脏计算模型陆续出现,并被成功应用于心脏功能研究与药物安全性评估[5-8]。
在早期虚拟生理心脏研究中,心脏一个生物秒的电生理活动往往需要数日甚至数月来仿真计算。随着显存技术的发展,这个时间缩短到数天。近年,有研究致力于提升虚拟生理心脏的计算速度。比如通过将三维心脏空间划分为矩形子区域来实现并行心脏模拟[9],使运算速度大大提升。另一项研究通过WebGL将高性能心脏模拟扩展到普通计算机上[10],甚至有GPU的手机也可以模拟三维心室的电动态。一些研究试图通过自适应时间步长来提高运行速度[11,12],结果表明,固定时间步长比自适应时间步长方法具有更好的效率[11]。
但这些研究仅能达到“准实时运算”,离真正意义上的“实时运算”,即仿真时间与生物时间比达到1:1,还有难以逾越的距离,更不用说仿真精度的提升带来的运算量爆炸式增长。高计算复杂度带来的海量运算,使得虚拟生理心脏模型难以实现实时计算,阻碍其大规模应用。
2
实时计算
为了在更高分辨率、更高精度和更大规模的心脏模型上实现实时仿真,智源研究院开发了具有精细细胞电生理与解剖结构的人心室模型。该模型包含了19种细胞生理状态变量和70多个公式,能够实现复杂的心脏电生理与病理仿真,为临床与医药工业应用提供丰富的场景。
为实现实时计算,智源对模型底层计算进行了深度优化。针对心脏仿真中计算强度大和I/O密集等瓶颈问题,智源充分结合A100平台的硬件特点,设计了多种优化策略,如量化和循环展开。这些措施有效降低了计算复杂度和I/O,使得在更大规模和更高复杂度的心脏模型上实现了180倍的速度提升。
最终,智源虚拟心脏仿真系统实现了对心脏电生理功能的实时仿真,达到生物时间与计算时间比为1:0.84。这一成果不仅提升了心脏仿真系统的性能,还为更广泛的医学研究和临床应用提供了强有力的支持,标志着心脏仿真技术的又一重大里程碑进展。
图注2:实时心脏计算概览图
2.1 技术路线
在GPU的架构设计中,顺序访问内存(如连续的数据访问)相较于随机访问具有更高的性能。此外,在执行顺序访问时,通常会采用预取技术提前加载数据,以进一步提高访问效率。
同时,在虚拟心脏模型中,大约有2/3的物理空间位置是空余腔体空间,有效心肌组织仅占1/3的物理空间。心脏仿真的主要计算和I/O操作都集中在对有效心肌组织中的每一个单细胞中的离子通道和细胞膜电位进行时间上的精细更新,同时考虑邻近细胞的电耦合影响。
基于GPU访存特点和心脏解剖结构的特殊性,我们设计了适合稀疏数据的数据结构。利用顺序访存提升I/O速度,确保并行线程仅处理有效细胞,从而最大限度地提高GPU内存的利用率。通过这种创新的结构,显著优化了计算性能,使得心脏仿真能够在IO访存上达到最优效果。
图注3:心脏模型有效数据在GPU内存上的排布
在计算层面,采用量化策略,有效简化模型中的对数和指数等复杂计算,从而显著降低了计算复杂度。
此外,为进一步减少I/O操作次数,采用循环展开策略,实现在一次读取中进行多次计算,大大降低I/O,显著提升SM核心的计算利用率。
基于A100平台,我们设计了高效的P2P通讯方式,利用GPU直连实现在节点内快速的数据交换,确保数据传输的低延迟与高带宽。在节点之间,采用RDMA(远程直接内存访问),进一步增强跨节点数据传输的效率,充分发挥硬件平台的并行计算与通讯能力。
图注4:技术路线图
2.2 仿真结果
我们测试了不同优化策略对仿真2生物秒心脏功能所用计算时间的影响,结果如下图所示。对2生物秒心脏功能的模拟,基准模型在未优化的情况下A100单卡需要计算时间为304.25秒。在采用分布式、量化、循环展开策略后,其所用时间分别是9.75、3.93、1.68秒。其中采用循环展开后,计算时间达到2秒内,达到计算时间/生物比小于1,实现实时/超实时计算的要求。其中,分布式计算对于系统仿真速度影响最大,达到了32倍提速。量化策略和循环展开策略分别将仿真速度提升了2.48和2.34倍。在同时采用分布式、量化、循环展开策略的情况下,系统仿真速度整体提升了181倍。
图注5:不同优化策略的计算时间
图注6:不同优化策略的速度提升
2.2.1 拓展曲线
图注7:不同优化策略的拓展曲线
如图7扩展曲线所示,随着GPU卡数的增加,基准模型和优化后的模型仿真时间都在减少。基准模型在增加到48卡后,计算时间不再减小。此时的生物:计算时间比为1:5。再采用量化和循环展开策略后,32张卡即可实现实时计算,生物:计算时间比达到1:0.84。
2.2.2 主要GPU指标
图注8:不同优化策略的计算密度和计算强度
图注9:不同优化策略的内存和SM利用率
通过GPU指标可以看出(图8,图9)量化策略通过提升IO同时降低计算的方式提高整体计算性能;循环展开通过大幅度降低I/O同时提高计算密度的方式提高计算性能。
2.2.3 计算精度
我们统计了加速前与加速后的结果误差,仿真的膜电位V的时程差别<2 ms (0.6%),模电位平均误差为0.72mV(0.4%),均满足生理准确度要求。优化前后主要离子通道的仿真曲线吻合(如图10所示)。
图注10:仿真前后细胞主要离子通道电流与胞内离子浓度在一心律节拍间的变化
3
总结
智源研究院从心脏模型的解剖结构、心肌细胞电生理的计算特点及计算系统的硬件架构出发,设计了心脏仿真系统的数据结构和优化策略,以提高计算效率。我们采用先进的并行处理方法,充分利用现代GPU设备的强大计算能力,优化数据传输和通讯方式,以减少延迟并提高数据吞吐量。通过这些策略,不仅提升了仿真系统的计算速度,还保证了在可接受误差范围内的计算精度,最终成功实现了心脏仿真的实时计算目标,达到超实时计算结果。这一成果为进一步研究心律失常产生的离子通道与分子机制等关键医学问题,也为手术规划如房颤射频消融方案等临床应用,以及新药研发与其心脏安全性筛选奠定了坚实基础,同时也为其它超大复杂物理系统的实时仿真提供坚实基础。
参考文献
1. Hodgkin A L, Huxley A F. A quantitative description of membrane current and its application to conduction and excitation in nerve[J]. The Journal of Physiology, 1952, 117(4): 500-544.
2. Noble D. Cardiac action and pacemaker potentials based on the Hodgkin-Huxley equations.[J]. Nature, 1960, 188(4749): 495-497.
3. Crampin E J , Halstead M , Hunter P , et al. Computational physiology and the physiome project[J]. Experimental Physiology, 2004, 89.
4. Smaill B, Hunter P. Structure and function of the diastolic heart: material properties of passive myocardium[M]//Theory of heart. Springer New York, 1991: 1-29.
5. Alday E A P, Colman M A, Langley P, Zhang H. Novel non-invasive algorithm to identify the origins of re-entry and ectopic foci in the atria from 64-lead ECGs: A computational study[J]. PLoS Computational Biolog. 2017, 13(3): e1005270.
6. Boyett M R, Li J, Inada S, et al. Imaging the heart: computer three-dimensional anatomical models of the heart[J]. Journal of Electrocardiology, 2005.
7. Nordsletten D, Niederer S, Nash M P, et al. Coupling multi-physics models to cardiac mechanics[J]. Progress in Biophysics & Molecular Biology, 2011, 104(1): 77-88.
8. Colman M A, Aslanidi O V, Kharche S, et al. Pro‐arrhythmogenic effects of atrial fibrillation‐induced electrical remodelling: insights from the three‐dimensional virtual human atria[J]. The Journal of physiology, 2013, 591(17): 4249-4272.
9. Wang W, Xu L, Cavazos J, Huang HH, Kay M. Fast acceleration of 2D wave propagation simulations using modern computational accelerators. PLoS One. 2014;9(1):e86484.
10. Kaboudian A, Cherry EM, Fenton FH. Real-time interactive simulations of large-scale systems on personal computers and cell phones: Toward patient-specific heart modeling and other applications. Sci Adv. 2019;5(3):eaav6019.
11. Garcia-Molla VM, Liberos A, Vidal A, Guillem MS, Millet J, Gonzalez A, et al. Adaptive step ODE algorithms for the 3D simulation of electric heart activity with graphics processing units. Comput Biol Med. 2014;44:15-26.
12. Sachetto Oliveira R, Martins Rocha B, Burgarelli D, Meira W, Jr., Constantinides C, Weber Dos Santos R. Performance evaluation of GPU parallelization, space-time adaptive algorithms, and their combination for simulating cardiac electrophysiology. Int J Numer Method Biomed Eng. 2018;34(2).