第十一届全国计算数学年会 #离子通道模拟 ####卢本卓 ####科学院计算数学与科学工程计算研究所, ####国家数学与交叉科学中心,北京 ####西安,07/21/2017
#内容 ## 生物背景 ## 电扩散PNP模型及其推广 ## 生物分子表面/立体网格生成 ## 有限元模拟 ## 应用举例 --纳米管的电流特性及第三代基因测序技术 --钾离子通道的选择性通透
#1. Ion channel - membrane system ![](http://data.xyzgate.com/10ab27051c52658aeec204a2609263ff.jpeg)
=!= A real view of protein molecule: <iframe frameborder="no" border="0" marginwidth="0" marginheight="0" width=100% height=400px src="/scene?id= 59f83a764f022e388f25e817"></iframe>
Feature: selective permeation Different Ions Carry Different Signals through Different Channels ![](http://data.xyzgate.com/053df684d3059e1a5d38d5da9edbc956.png) Thanks Bob Eisenberg for the slide!
##All-atom MD simulation ![](http://data.xyzgate.com/c387b5b9037f19b17dd769b384210dcf.png)
#2. Continuum models **Poisson-Nernst-Planck** type of model: ```math \begin{cases} \frac{\partial c^i(r,t)}{\partial t}=\nabla \cdot {D^i(\nabla c^i + \beta c^i \nabla(q^i \phi ))},\ i=1,\cdots,K,\ \\ \nabla \cdot \varepsilon \nabla \phi(r,t) =-\rho(x) =-\displaystyle\sum q_m \delta(r-r_k)-\lambda \displaystyle\sum q^i c^i(r,t) \end{cases} ``` ![](http://data.xyzgate.com/391214b1ad3d1fcc7ad4308cb386000c.png) ![](http://data.xyzgate.com/6e6bc51318a8d619cd85c0a824f5e298.png) ```math \Rightarrow I=\int_s \displaystyle\sum_i^K D^i (\nabla c^i +\beta c^i \nabla(q^i \phi)) ```
##Some model extensions Derive from a complete free energy functional form: ```math \begin{aligned} F[c] = \int_{\Omega}\{\frac{1}{2}\rho(c)\phi(c) + \beta^{-1}\sum_{i=1}^{K}c_{i}[\log(\Lambda^{3}c_{i})-1] + \sum_{i=1}^{K}\mu_{i}^{ex}c_{i} \}dV + \int_{\Gamma_{N}}\frac{1}{2}\sigma\phi(c)dS -\int_{\Gamma_{D}}\frac{1}{2}\epsilon(c)\frac{\partial \phi(c)}{\partial n}\phi_{0}dS, \end{aligned} ``` *X. J. Liu, Y. Qiao, and B. Z. Lu. https://arxiv.org/abs/1702.01115v1 * - **Size effects** --> **SMPNP** (Lu B. *Biophysical J*. 2011 ) ```math \nabla \cdot \varepsilon(r) \nabla \phi(r)=-\rho^f(r)-\displaystyle\sum_i q^i c^i(r), \ \ r\in \Omega, ``` ```math -\nabla \cdot \{ {D^i(r) \nabla c^i(r)+ \frac{D^i(r)v^i(r)c^i(r)}{1-\displaystyle\sum_k a_k^3 c^k(r)} \displaystyle\sum_k a_k^3 \nabla c^k(r)+ \beta D^i(r)c^i(r)q^i \nabla \phi(r)} \}=0,\ \ r\in \Omega, \ i=1, \cdots, K, ``` - **When epsilon(c) --> Variable dielectric coefficient model VDPNP ** ```math -\nabla\cdot(\epsilon(c)\nabla\phi(c)) = \rho^{f} + \sum_{i=1}^{K}q_{i}c_{i}, \quad in \;\Omega, ``` ```math \frac{\partial c_{i}}{\partial t}=\nabla \cdot (D_{i}[\nabla c_{i}+\beta c_{i}\nabla(q_{i}\phi(c)- \frac{1}{2}\epsilon_{i}^{'}(c)\nabla\phi(c)\cdot \nabla\phi(c))]), \quad in \;\Omega_{s}, i=1, 2, \cdots, K. ``` - **Ion solvation energy --> BPNP ** (see later)
#3. Molecular mesh generation ![](http://data.xyzgate.com/1a43bce46d141a1b69b519a69952931e.png) ![](http://data.xyzgate.com/570963cd1ebf7dec005dd9d0e3c4b3d7.png)
## New surface meshing in our TMSmesh2 输入:PQR 文件,包含分子内所有原子的中心和半径。 输出:OFF 文件,包含生成的三角形表面网格。 处理的分子表面:高斯表面 ```math \{ \vec x \in R^3,\phi \left(\vec x \right) = c \}, ``` ```math \phi \left(\vec x \right) = \sum _{i=1}^{N}e^{-d(\Vert \vec x - \vec x_i \Vert ^2 -r_i^2)}, ``` Duncan, B. S.; Olson, A. J. *Biopolymers* **1993**, 33, 231–238. 算法分为两个步骤 1.第一步是自适应估计和空间划分过程。在这一步中,空间被自适应的划分为一些立方体盒子,在每个立方体盒子中用三线性表面来近似高斯表面。 2.第二步是把每个立方体盒子中的三线性表面剖分成沿 $$x, y, z$$ 方向单值的面片,然后进行三角化。 Tiantian Liu, Minxin Chen, Benzhuo Lu, Efficient and Qualified Mesh Generation for Gaussian Molecular Surface Using Adaptive Partition and Piecewise Polynomial Approximation, https://arxiv.org/pdf/1611.03474.pdf (submitted to SIAM Journal on Scientific Computing).
## n 阶多项式近似高斯表面 ```math \phi \left(\vec x \right) = \sum _{i=1}^{N}e^{-(\Vert \vec x - \vec x_i \Vert ^2 - r_i^2 )}=\sum _{i=1}^{N}e^{r_i^2}e^{-\left( x-x_i \right)^2}e^{-\left( y-y_i \right)^2}e^{-\left( z-z_i \right)^2}. ``` 在范围是 $$[a,b]\times [c,d]\times [e,f]$$ 的立方体盒子中,方程可以近似为: ```math P(x,y,z) = \sum_{i=1}^{N}e^{r_i^2}P_n(x,x_i,a,b)Q_n(y,y_i,c,d)R_n(z,z_i,e,f), ``` ```math P_n(x,x_i,a,b) = \sum_{j=0}^n \alpha_j(x_i,a,b)L_j(\frac{2x-(a+b)}{b-a}), ``` ```math Q_n(y,y_i,c,d) = \sum_{j=0}^n \beta_j(y_i,c,d)L_j(\frac{2y-(c+d)}{d-c}), ``` ```math R_n(z,z_i,e,f) = \sum_{j=0}^n \gamma_j(z_i,e,f)L_j(\frac{2z-(e+f)}{f-e}), ```
##三线性多项式近似高斯多项式 1.把一个可能存在交点的立方体盒子细分成 8 个子盒子,每个子盒子中对应的多项式依然可以用勒让德多项式来表示。 2.当子盒子越来越小时,系数张量中高阶勒让德多项式对应的系数越趋近于 0。 3.经过多次的细分和估计之后,在每个最终得到的立方体盒子里面,用三线性表面来近似多项式表面 $$\tilde P(x,y,z) = c$$。 ```math g(x,y,z) = a_0xyz+a_1xy+a_2xz+a_3yz+a_4x+a_5y+a_6z+a_7 = c. ```
##自适应估计和空间划分示例 分子ADP的网格生成过程中的立方体盒子 ![](http://data.xyzgate.com/00e2722de0c8d321191070bf32dc405d.png) ##三线性表面的三角化 ![](http://data.xyzgate.com/d0e970ce3d66ac66f4436e32fdde81f7.png)
## Example results 分子2JM0 ![](http://data.xyzgate.com/10f28ab4847eaddf859f250508dd2616.png) 30S ribosome ![](http://data.xyzgate.com/4b5096db60d0e6e8b5a1bf81752948ac.png)
##Volume mesh generation - Using Tetgen to generate volume mesh based on TMSmesh surface mesh - Sometime it needs to further reduce/smooth the surface mesh before volume mesh generation MX Chen, B Tu, and BZ Lu, Surface Triangular Mesh and Volume Tetrahedral Mesh Generations for Biomolecular Modelling, in book *Image-Based Geometric Modeling and Mesh Generation*, Publisher, Springer Berlin Heidelberg, 2013, XI, p85-106.
##Meshing for membrane-protein system - Difficulty: distinguish the tetrahedra in membrane region and pore region, which may be connected by holes or crevices. ![](http://data.xyzgate.com/b7e126d09612416c42f1a25eb95776f9.png) Walk-and-detect algorithm ![](http://data.xyzgate.com/569a76aa36eb5d7dc8162696e0995f78.png) TT Liu, SY Bai, B Tu, MX Chen, and BZ Lu, Membrane-Channel Protein System Mesh, Construction for Finite Element Simulations, *Mol. Based Math. Biol*. 2015; 3:128–139
## Connexin (Cx26) ![](http://data.xyzgate.com/6a5af8351ad734e438ebbcbfe7dad2a1.png) ![](http://data.xyzgate.com/562252ebdc0ef79d325b4099ea681eed.png) ![](http://data.xyzgate.com/a2218918d65bf73335d31cafaf75e0a9.png)
#4. FEM numerial algorithms - Parallel adaptive Finite element method (based on PHG) - Interface problem (jump dielectric coefficient) -- Conforming mesh (body-fitted mesh generation) - Singular charges treatment -- Potential decomposition to remove the singularities - Solving the coupled systems (PNP, smPNP …) -- Gummel iteration -- Stablization method (SUPG) - **Main issues**: -- Stable and Robust algorithm (sometimes the isuue caused by mesh quality) -- efficiency * Zhang LB, Numer. Math. Theor. Meth. Appl. 2, 65 (2009). Cheng IL et al. 2003 JH Chaudhry, J Comer, A Aksimentiev, LN Olson, Commun. Comput. Phys. 15: 93, 2014 Lu BZ, et al. 2007; 2008; 2009; 2010; 2011; 2013,2015 *
#5. Some FEM simulations and applications **Current-voltage characteristics** ## DNA-nanopore/channel sequencing ### Cylindrical pore ![](http://data.xyzgate.com/2386ff702b74793988ac1e401bd75d46.png) ![](http://data.xyzgate.com/a5961765fa329d820cb61cd0263609de.png) A dsDNA in a 4-nm-diameter nanopore. The SiN membrane is shown in gray. (submitted)
### Influence of membrane thickness ![](http://data.xyzgate.com/e40ced820f46034d6b0d9f09aeccd477.png) ![](http://data.xyzgate.com/f19726fd9222dfa753b436487344a805.png) Left: Dependence of average experimental Io (without DNA, black circles) and the most probable DNA current amplitude Ip (with DNA, red triangles) on h. The black dashed line is a fit using equation (1) to the average Io data from the combined data of 20 pores, which yields an effective pore thickness. Right: computational results.
##Conic pore works better! ###--Influence of conic pore diameter and shape ![](http://data.xyzgate.com/df9ac5604fd11f30e49b74b27abbd5d3.png) ![](http://data.xyzgate.com/5135b2ca0e2899606c4c17b035f189ab.png) Open-pore currents Is and current amplitudes ∆Is with 20-pb dsDNA in the conical nanopores. (2r = 4 nm)
## K+ channel selectivity - Using Born-energy-PNP model - KcsA structure: PDB code 1BL8. *MacKinnon et al, Science. 280,(1998) 69.* ![](http://data.xyzgate.com/7d8945b8a2fa39d50f58ecb7652b400b.png) ![](http://data.xyzgate.com/d6cb8b9a0a743820dfc2d69a296a84cb.png)
### Considering Born energy - Born solvation energy ```math u_{ex}:\ G_{born}=\frac{q_i^2}{a_i}(\frac{1}{\varepsilon(r)}-\frac{1}{\varepsilon_o}) ``` ![](http://data.xyzgate.com/8bf93fa3ecb96a2d2127972c16ac7fc4.png) $$\to$$ **BPNP** ```math \nabla \cdot \varepsilon(r) \nabla \phi(r)=-\rho^f(r)-\displaystyle\sum_i q^i c^i(r), \ \ r\in \Omega, ``` ```math \nabla\cdot[D^i(r)(\nabla c^i(r)+ \beta c^i(r) \nabla(q_i \phi(r)+ \alpha \frac{q_i^2}{2a_i}( \frac{1}{\varepsilon(r)-\frac{1}{\varepsilon_o}}))=0, \ \ r\in \Omega_s, \ \ i=1, \cdots, K, ``` *X. J. Liu, Y. Qiao, and B. Z. Lu. https://arxiv.org/abs/1702.01115v1 *
### 网格 ![](http://data.xyzgate.com/024a4c13552ac663d7a55feced752edd.png)
### 介电系数 ![](http://data.xyzgate.com/f7c6d1364160668ed02b9d725019f081.png)
### K+选择性及内向整流 Conditions: Mixed electrolyte: Na+, K+, Cl- $$c_\text{bulk}$$ for Na+ and K+ = 0.1M, $$c_\text{bulk}$$ for Cl- = 0.2M, radii:$$a_{Na} = 1.0A, a_K = 1.50A and a_{Cl} = 2.0A$$. Ion distributions under a fixed membrane voltage φ0 = 0.20V. ![](http://data.xyzgate.com/8e79153e7986e0e204125c21e32b9054.png)
I-V curves ![](http://data.xyzgate.com/f98fd08520fa12d984106a7b4f518675.png) (*paper to be submitted*)
#6. Web-server www.xyzgate.com
# Take-home message 离子通道 - 选择性 - 模型 - 算法 - 网格生成
# Acknowledgments 计算数学所: 白石阳,谢妍,乔瑜,刘田田,刘雪娇,许竞劼 张林波 研究员 陈旻昕 (苏州大学) 涂斌 (国家纳米中心) Robert Eisenberg (Rush Univ., US) 经费:CAS, NSFC, NCMIS等 #Thanks !
##Some model extensions Free energy ```math F = \int [\frac{1}{2} \varepsilon \phi(x) \rho(x)dx +k_B T \displaystyle\sum_i c_i (\ln c_i \Lambda^{3} -1) + \displaystyle\sum_i c_i u_i^{ex}]dx ``` - **Size effects** --> **SMPNP** (Andelman et al, PRL. 79:435, 1997; Lu B. *Biophysical J*. 2011 ) ```math \nabla \cdot \varepsilon(r) \nabla \phi(r)=-\rho^f(r)-\displaystyle\sum_i q^i c^i(r), \ \ r\in \Omega, ``` ```math -\nabla \cdot \{ {D^i(r) \nabla c^i(r)+ \frac{D^i(r)v^i(r)c^i(r)}{1-\displaystyle\sum_k a_k^3 c^k(r)} \displaystyle\sum_k a_k^3 \nabla c^k(r)+ \beta D^i(r)c^i(r)q^i \nabla \phi(r)} \}=0,\ \ r\in \Omega, \ i=1, \cdots, K, ```
- Ion solvation energy --> BPNP
When $$\epsilon(c)$$ --> - **Variable dielectric coefficient model VDPNP ** Note: a complete free energy functional is presented: ```math \begin{aligned} F[c] = \int_{\Omega}\frac{1}{2}\rho(c)\phi(c) dV + \int_{\Gamma_{N}}\frac{1}{2}\sigma\phi(c)dS -\int_{\Gamma_{D}}\frac{1}{2}\epsilon(c)\frac{\partial \phi(c)}{\partial n}\phi_{0}dS + \beta^{-1}\sum_{i=1}^{K}\int_{\Omega}c_{i}[\log(\Lambda^{3}c_{i})-1]dV - \sum_{i=1}^{K}\int_{\Omega}\mu_{i}c_{i}dV, \end{aligned} ``` Now we get a set of generalized self-consistent PNP equations with concentration-dependent dielectric (VDPNP): ```math -\nabla\cdot(\epsilon(c)\nabla\phi(c)) = \rho^{f} + \sum_{i=1}^{K}q_{i}c_{i}, \quad in \;\Omega, ``` ```math \frac{\partial c_{i}}{\partial t}=\nabla \cdot (D_{i}[\nabla c_{i}+\beta c_{i}\nabla(q_{i}\phi(c)- \frac{1}{2}\epsilon_{i}^{'}(c)\nabla\phi(c)\cdot \nabla\phi(c))]), \quad in \;\Omega_{s}, i=1, 2, \cdots, K. ``` ```math \epsilon(c)\frac{\partial\phi}{\partial n} = \sigma \quad on \;\Gamma_{N}, ``` ```math \phi = \phi_{0} \quad on \;\Gamma_{D}, ``` ```math c_{i} = c_{i}^{b} \quad on \;\Gamma_{D}, ``` ```math J_{i}\cdot n = 0 \quad on \;\Gamma_{m}. ``` If the dielectric coefficient does not depend on local ionic concentrations, VDPNP reduces to the traditional PNP equations. *X. J. Liu, Y. Qiao, and B. Z. Lu. https://arxiv.org/abs/1702.01115v1 *
## Diffusion coefficient inside the the channel ```math D(r)=\begin{cases} D_{bulk}, \ \ r\in bulk\ region, \\ D_{chan}+(D_{chan}-D_{bulk})f(r), \ \ r\in buffering\ region \\ D_{chan},\ \ \ r\in channel\ region, \end{cases} ``` ```math f(r)=f(z)=n(\frac{z-z_{chan}}{z_{bulk}-z_{chan}})^{n+1}-(n+1)(\frac{z-z_{chan}}{z_{bulk}-z_{chan}})^n ``` Hwang, H.; Schatz, G. C.; Ratner, M. A. *J. Phys. Chem.* *A **2007**, 111,* 12506–12512 where D_{Cl} = 0.203Å^2/ps, D_K = 0.196Å^2/ps, n = 9. For the bottom boundary, $$z_{chan} = 11$$ and $$z_{bulk} = 9$$. For the top boundary, $$z_{chan} = 32 and z_{bulk} = 34$$. B Tu, BZ Lu, et al. *J. Comput. Chem.* 2013, 34: 2065
/
Go to slide: