Parallel Numerical Solution of Intracellular Calcium Dynamics - PowerPoint PPT Presentation

Parallel Numerical Solution of Intracellular Calcium Dynamics Chamakuri Nagaiah 1 , Gerald Warnecke 1 udiger 2 , Martin Falcke 2 Sten R¨ 1 Institute for Analysis and Numerics Otto-von-Guericke University, Magdeburg 2 Department of Theoritical Physics Hahn-Meitner Institute, Berlin July 3-7 2006 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 1 / 26

Outline Introduction 1 Governing Equations in 2D 2 Grid Adaptivity 3 FEM Discretization 4 Numerical Results 5 Conclusions and Future Work 6 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 2 / 26

Introduction The endoplasmic reticulum (ER) has a high calcium concentration Channels: ER → cytosol, pumps: cytosol → ER Ca 2 + is released by transient openings of channels Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 3 / 26

Introduction channel opens, releases Ca 2 + from the ER into the cytosol Ca 2 + diffuses to neighboring channels increase of Ca 2 + favors opening: amplification very high Ca 2 + decreases opening probability: inhibition Ca 2 + is pumped back from the cytosol into the ER Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 4 / 26

Structure of Cluster and Channels 20,000-1000,000 nm Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 5 / 26

Experimental Result < experimentmovie . avi > Puffs and waves in the stochastic regime (I. Parker, UC Irvine) Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 6 / 26

Deterministic Equations in 2D c 2 ∂ c � = D c ∆ c + ( P l + P c ( r ))( E − c ) − P p d + c 2 − H i ( c , b i ) K 2 ∂ t i c 2 ∂ E � � � = D E ∆ E − γ ( P l + P c ( r ))( E − c ) − P p − K j ( c , b E , j ) K 2 d + c 2 ∂ t j ∂ b i = D b , i ∆ b i + H i ( c , b i ) , i = 0 , n − 1 ∂ t ∂ b E , j = D E , j ∆ b E , j + K j ( E , b E , j ) , j = 0 , m − 1 . ∂ t where H i = k + b , i ( B i − b i ) c − k − b , i b i E , j ( G i − b E , j ) E − k − K j = k + E , j b E , j B.C’s: no flux at the boundaries [Thul 04; Falcke 03,04]. Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 7 / 26

Determination of P c ( r ) Each cluster is given by a fixed position � X i and its radius � R i = R S N open , i N open , i is the number of open channels in cluster i . This number is determined by channel dynamics. � � � r i − � � � P ch if X i � < R i for a cluster i � � P c ( � r i ) = 0 otherwise Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 8 / 26

Zienkiewicz-Zhu Error Indicator Let Gu h ∈ V h be the � ., . � h -projection of ∇ u h onto V h , calculated by | T | � Gu h ( x i ) = | w x |∇ u h | T ( x i ) T ⊂ w x Error estimator η Z , T := � Gu h − ∇ u h � L 2 ( T ) and 1 / 2     � η 2 η Z := Z , T  T ∈T h  O. C. Zienkiewicz, J. Z. Zhu. A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. Num. Meth. Eng . 24 (1987) 337-357 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 9 / 26

Zienkiewicz-Zhu Error Indicator Let λ ( T ) ∈ N 0 be the refinement level of triangle T ∈ T h , λ max ∈ N 0 be a given maximum refinement level ϕ 1 , . . . , ϕ λ max given real numbers satisfying 0 ≤ ϕ 1 ≤ . . . ≤ ϕ λ max. Triangle T is marked for refinement if η Z ∈ [ ϕ i , ϕ i + 1 ] and λ ( T ) < i , i = 0 , . . . , λ max. Adaption parameters are λ max = 6 and ϕ 1 = 0 . 0001 , ϕ 2 = 0 . 0002 , ϕ 3 = 0 . 0004 , ϕ 4 = 0 . 0008 , ϕ 5 = 0 . 0016 , ϕ 6 = 0 . 0032. Programm package UG (Unstructured Grid) developed by group of G. Wittum and P . Bastian at University of Heidelberg. P . Bastian, K. Birken, S. Lang, K. Johannsen, N. Neuß, H. Rentz-Reichert and C. Wieners. UG: A flexible software toolbox for solving partial differential equations. Computing and Visualization in Science, 1 (1997) 27–40 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 10 / 26

Adaptive Mesh for Single Cluster level 0 (initial mesh) level 1 level 6 nodes = 2378 nodes = 2433 nodes = 2766 elements = 4566 elements = 4676 elements = 5342 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 11 / 26

Adaptive Mesh for 100 Clusters level 0 (initial mesh) level 1 level 6 nodes = 3503 nodes = 4964 nodes = 19367 elements = 6776 elements = 9698 elements = 38204 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 12 / 26

Typical Problem ∂ u ( x , t ) − ∇ · ( a ( x ) ∇ u ( x , t )) + s ( u ( x , t )) = f ( x , t ) in Ω × ( 0 , T ] ∂ t ∂ u ( x , t ) = 0 on ∂ Ω × ( 0 , T ] . ∂η Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 13 / 26

Spatial Discretization Weak formulation find u ∈ V = H 1 (Ω) such that � ∂ u ∂ t , v � + � a ( x ) ∇ u , ∇ v � + � s ( u ) , v � = � f , v � for all v ∈ V � ∂ u h ∂ t , v h � + � a ( x ) ∇ u h , ∇ v h � + � s ( u h ) , v h � = � f , v h � for all v h ∈ V h Approximate the solution u h ∈ V h using the basis functions u h ( t , x ) = � N i = 1 u i ( t ) φ i ( x ) Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 14 / 26

Spatial Discretization M˙ u + Au + S = F (1) where M - mass matrix, A - stiffness matrix. The matrices are defined as follows, M = � φ i , φ j � , A = � a ( x ) ∇ φ i , ∇ φ j � S = � s ( � N i = 1 u i ( t ) φ i ( x )) , φ j � F = � f , φ j � . Approximate the term S using quadrature rule S = � φ i , φ j � s ( u i ) = MR . Mass lumping, inverting the lumped mass matrix u = − M − 1 Au − R + M − 1 F ˙ (2) Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 15 / 26

Time Discretization We considered the ODE problem ∂ u u ( t 0 ) = u 0 . ∂ t = F ( t , u ) , (3) The i -th time step of a W-method of order p with embedding of order ˆ p � = p has the form j − 1 j − 1 � � t i + τ i a j , u i + τ i � � ( I − τ i γ J ) k j = + c lj k l , j = 1 , . . . , s , (4) F b lj k l l = 1 l = 1 s u i + τ i � u i + 1 = d l k l , (5) l = 1 s u i + τ i ˆ u i + 1 � ˆ = d l k l . (6) l = 1 B. A. Schmitt and R. Weiner. Matrix-free W-methods using a multiple Arnoldi iteration. Appl. Num. Math. 18 (1995) 307-320 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 16 / 26

Time Discretization We use a W-method with s = 3 stages given by the coefficients √ γ = 1 − 1 2 , a 1 = 0 , a 2 = 1 , a 3 = 1 , b 12 = 1 , b 13 = 1 , b 23 = 0 , 2 √ √ c 12 = − 2 − 2 , c 13 = − 1 , c 23 = − 1 + 2 , √ d 2 = 1 2 − 1 d 3 = 1 d 1 = 1 , 2 , 2 , 2 √ √ √ d 1 = 9 10 − 1 d 2 = 9 20 − 11 d 3 = 11 20 + 1 ˆ ˆ ˆ 2 , 2 , 2 . (7) 20 20 20 This method was used by Schmitt and Weiner for the construction of a Krylov-W-method. Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 17 / 26

Time Discretization A new time step τ new is computed by  β max τ i , τ > β max τ i , ¯ 1 � TOL t � p + 1 ˆ  τ := βτ i β min τ i , τ < β min τ i , ¯ , τ new := ¯ ǫ τ, ¯ otherwise .  where β > 0 is safety factor. Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 18 / 26

Linear Solver Solve s linear systems of the general form ( I − τ i γ J ) k j = b j , j = 1 , . . . , s . where k j are the unknown vectors and A := I − τ i γ J is the same for all stages. Iterative solvers are appropriate BiCGSTAB method with ILU preconditioning. Tolerance for the linear solver is TOL LS = α LS TOL t /τ i H.A. van der Vorst. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 13 (1994) 631-644 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 19 / 26

Numerical Results < Calcium 100 clustersStoc . avi > Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 20 / 26

Time Step Rejections Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 21 / 26

Numerical Result with One Opening Channel Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 22 / 26

Domain Decomposition using RIB algorithm B. Hendrickson, R. Leland. The CHACO user’s guide 1.0. Technical Report SAND93-1301, Sandia National Laboratory, 1993. Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 23 / 26

Domain Decomposition using RIB algorithm nodes = 33,322 elements = 66,370 unknowns = 133,288 nodes = 32,417 elements = 64,560 unknowns = 129,668 Chamakuri Nagaiah (IAN) Otto-von-Guericke Universit¨ at Magdeburg DD17, Austria, July 3-7 2006 24 / 26