在复杂几何中的物质运输以保证功能是一个复杂的问题。通过等几何(IGA)模拟物质在生物神经网络中的运输是有意义的。
文章使用反应-扩散-运输(reaction-diffusion-transport)方程表示物质传输过程,用截断分层三次B样条(THB-Spline3D)表示神经的几何形状。并通过Navier-Stokes(NS)方程获得物质传输的速度场,最后使用IGA求解器模拟了单个轴突、二分叉轴突与三分叉轴突的物质传输过程。
1. 总览
文中提出的对神经网络物质传输的建模分析流程图如下所示
分为两部分:几个建模部分与仿真分析部分
在建模部分,采用基于骨架的扫掠方法来生成轴突的六面体控制网格,然后在控制网格上构建THB-Spline3D用来表示模拟几何体。
在仿真分析部分,通过将一维的辅助动力传输模型推广到三维,得到IGA求解器进行仿真求解。
2. 控制网格生成
为表达轴突中复杂分支的几何形状,使用中心线上的折点以及方向和相关直径等数据可视化模型。如图A所示,将每一段神经网络的中心线当做骨架,通过沿骨架扫掠横截面的四边形网格来生成每个分段的六面体控制网格。
生成网格的目的是指导样条线的生成,以便更准确、更平滑地表示轴突网络几何体。但由于一对一的扫掠要求源曲面与目标曲面具有相同的拓扑,而在轴突的分叉处并不满足这一特性,于是需要在每个分叉处加入一个半平面,如图B、C所示。
生成的网格要在每个分支点附近达到G1连续性,需要满足一下两点需求:1. 控制网格中两个分支共用的顶点须沿着轴向与相邻的两个折点共线;2. 三个分支共享的边界顶点必须与其所有相邻的边界顶点共面。使用上述的基于骨架的扫掠方式生成的几个简单模型控制网格如下:
接下来的工作就是使用得到的控制网格构造THB-Spline3D,以最终表示模型。最后使用得到的模型进行物质传输模拟。
2.1 THB-Spline3D综述
截断分层三次B样条(THB-Spline3D)是通过将四边形网格上的双三次混合函数推广到六面体网格上的三次三次混合函数得到的。在三次混合函数的基础上,进一步完成了三次混合函数的层次结构和截断机制,实现了高度局部化的细化。
在六面体网格中,我们定义一条边如果被周围4个六面体共用,则该边为非奇异边,否则为奇异边。而相对应的,若组成六面体的所有边都为非奇异边则该六面体为规则的,否则是不规则的。将六面体网格中的单元(六面体)分为3类:边界单元、内部规则单元、与内部不规则单元。
下图A描述了局部索引为1~8的六面体单元;图B描述了与该六面体单元对应的64个Bézier点,其中红色表示体点,黄色表示面点,绿色表示边缘点,蓝色表示角点。
在不规则的单元上应用Bézier提取矩阵 C \mathbf{C} C可以将局部的样条控制网格转化为Bézier控制网格。对于给定的六面体单元 Π e \Pi_e Πe,64个Bézier点可以通过其局部控制网格中的顶点 P \mathbf{P} P计算得到:
体点 R b o d y e \mathbf{R}_{body}^e Rbodye:通过8个顶点 P b o d y e \mathbf{P}_{body}^e Pbodye的凸组合, R b o d y e = C b o d y e P b o d y e \mathbf{R}_{body}^e=\mathbf{C}_{body}^e\mathbf{P}_{body}^e Rbodye=CbodyePbodye。 C b o d y \mathbf{C}_{body} Cbody将8个角点转变成8个Bézier体点。
非体点 R n o n − b o d y e \mathbf{R}_{non-body}^e Rnon−bodye:得到所有六面体的Bézier体点后,计算所有体点与非体点之间的欧氏距离。对于每个非体点都可以找到距离最近的体点,得到集合 Ω ( e , i ) \Omega(e,i) Ω(e,i)。之后可以对这些最近的体点求平均得到单个六面体的非体点 R n o n − b o d y e \mathbf{R}_{non-body}^e Rnon−bodye
R n o n − b o d y e , i = 1 n i ∑ ( k , j ) ∈ Ω ( e , i ) R b o d y k , j (1) R_{n o n-b o d y}^{e, i}=\frac{1}{n_{i}} \sum_{(k, j) \in \Omega(e, i)} R_{b o d y}^{k, j} \tag{1} Rnon−bodye,i=ni1(k,j)∈Ω(e,i)∑Rbodyk,j(1)
其中 n i n_i ni为共享该非体点 R n o n − b o d y e \mathbf{R}_{non-body}^e Rnon−bodye的六面体个数, ( k , j ) (k,j) (k,j)表示六面体 Π k \Pi_k Πk的第 j j j个非体点。
之后可以通过 R e = C P \mathbf{R}^e=\mathbf{CP} Re=CP在局部样条控制网格中导出64个贝塞尔点的凸组合, Π e \Pi_e Πe上的混合函数定义为 N e = C T b N^e=\mathbf{C}^Tb Ne=CTb,其中
b = [ g 0 ( ξ ) g 0 ( η ) g 0 ( ζ ) , g 1 ( ξ ) g 0 ( η ) g 0 ( ζ ) , … , g 3 ( ξ ) g 0 ( η ) g 0 ( ζ ) , g 0 ( ξ ) g 1 ( η ) g 0 ( ζ ) , g 1 ( ξ ) g 1 ( η ) g 0 ( ζ ) … , g 3 ( ξ ) g 3 ( η ) g 0 ( ζ ) , g 0 ( ξ ) g 0 ( η ) g 1 ( ζ ) , … , g 1 ( ξ ) g 0 ( η ) g 1 ( ζ ) , … , g 3 ( ξ ) g 3 ( η ) g 3 ( ζ ) ] T (2) \begin{aligned} b=&\left[g_{0}(\xi) g_{0}(\eta) g_{0}(\zeta), g_{1}(\xi) g_{0}(\eta) g_{0}(\zeta), \ldots, g_{3}(\xi) g_{0}(\eta) g_{0}(\zeta), g_{0}(\xi) g_{1}(\eta) g_{0}(\zeta), g_{1}(\xi) g_{1}(\eta) g_{0}(\zeta)\right.\\ &\left.\ldots, g_{3}(\xi) g_{3}(\eta) g_{0}(\zeta), g_{0}(\xi) g_{0}(\eta) g_{1}(\zeta), \ldots, g_{1}(\xi) g_{0}(\eta) g_{1}(\zeta), \ldots, g_{3}(\xi) g_{3}(\eta) g_{3}(\zeta)\right]^{T} \end{aligned} \tag{2} b=[g0(ξ)g0(η)g0(ζ),g1(ξ)g0(η)g0(ζ),…,g3(ξ)g0(η)g0(ζ),g0(ξ)g1(η)g0(ζ),g1(ξ)g1(η)g0(ζ)…,g3(ξ)g3(η)g0(ζ),g0(ξ)g0(η)g1(ζ),…,g1(ξ)g0(η)g1(ζ),…,g3(ξ)g3