High performance computing in multiscale modeling cardiac contraction: Bridging proteins to cells to whole heart

1. High performance computing in multiscale modeling cardiac contraction: Bridging proteins to cells to whole heart J. Jeremy Rice Computational Biology Center Thomas J. Watson Research Center johnrice@us.ibm.com

2. Motivation • Deaths due to cardiovascular diseases remain the largest contributor to premature mortality in most developed societies. • The cost of heart disease and stroke in the US is projected to be more than $448 billion in 2008. • Amiodarone is the most effective antiarrhythmic drug - but uncertain action and serindipitous discovery. • Predictive mathematical models are ideal tools to address the complex nature cardiovascular diseases. • Efforts to "customize" cardiac models are pushing towards role as enablers of personalized medicine.

3. Building Multiscale Heart Models Heart will allow better therapies for heart disease.. ...but will require bridging between organ level and molecular level Organ level Cell level Molecular level Sarcomere contracts by cyclical interactions of myosin on thick filament (red) and actin in thin filament (green). In each cell of heart, a lattice of sarcomeres produce contraction on every heart beat. Reconstruction of whole heart by Peter Hunter, U. of Auckland

4. Cyclical actin and myosin interactions converts energy in ATP to force/motion From the website of the Michael Geeves' Laboratory: http://www.kent.ac.uk/bio/geeves/Research/home.htm

5. Banded Structure of Sarcomere

6. Ca Ca Ca Building a multiscale model on actin-myosin in muscle B. Sarcomere component proteins A. Sarcomere structure z-line z-line thin filament troponin complex … … actin TnC TnT TnI … … … … tropomyosin myosin head thick filament myosin neck z-line intertwined myosin tail regions Single thick and thin filament in half sarcomere

7. Thin filament is a two-stranded helix of actin monomers "Pseudo-repeat" 37 nm "Pseudo-repeat" = 13 units True repeat = 26 units 5.54 nm 2.77 nm From http://www.kent.ac.uk/bio/geeves/Research/home.htm

8. Position of myosin heads on thick filament • Pairs of heads emanate 180 degree apart in radial direction at each step • Radial direction of heads rotate ~60 degrees at next step in axial direction (distance = ~14.3 nm) • a "pseudo-repeat" happens on the 3th steps as heads will be radiate in same radial direction (distance = ~43 nm) axial direction radial direction

9. A. Spatial relationships of actin and myosin yj+1 yj 37 nm C. Markov state model of crossbridging … … 43 nm xj+1 xj B. Mechanical representation ka= 1743 pN/nm yj+2 yj yj-1 yj+2 … … -7 nm kxb= 1 pN/nm … … xj+2 xj+1 xj xj-1 km= 2020 pN/nm

10. Representing compliances in myofilament Daniel, Trimble & Chase, 1998 K.X = A

11. Towards a more complete computational model

12. Cardiac Muscle Physiology

13. Cardiac Muscle Physiology

14. Ca Ca Representing nearest-neighbor interactions kon[Ca] 0N 1N koff gnknp_1 g-nkpn _1 g-nkpn_0 gnknp_0 k’on[Ca] 1P 0P k’off g > 1 for neighbor-to-neighbor cooperativity n = number of activated neighbors in 0P or 1P state (0, 1 or 2)

15. Compute exponent on g … … n = 0 … … n = 1 { or … … … … n = 2

16. 1N 0N Ca Ca Ca 1N 0N Ca 1P 0P Ca Regulation Ca 1P 0P Ca 1PreF 0PreF P 0F 1F PreF F Ca regulation and XB Cycle XB Cycle Combining computational models

17. Ca Ca Ca Ca Ca Ca Ca Must Use Monte Carlo Methods Relaxed Activated – Generates Force

18. Sample Results - Variation of Gamma Ca Ca n = 0 g Exponent on n = 1 n = 2

19. Some Results - F-Ca Relations Unpublished data from Dobesh et al., 2001, rat skinned fiber, SL = 2.15 mm

20. Ca Ca Representing nearest-neighbor interactions kon[Ca] 0N 1N koff gnknp_1 g-nkpn _1 g-nkpn_0 gnknp_0 k’on[Ca] 1P 0P k’off g > 1 for neighbor-to-neighbor cooperativity n = number of activated neighbors in 0P or 1P state (0, 1 or 2)

21. Compute exponent on g … … n = 0 … … n = 1 { or … … … … n = 2

22. 1N 0N Ca Ca Ca 1N 0N Ca 1P 0P Ca Regulation Ca 1P 0P Ca 1PreF 0PreF P 0F 1F PreF F Ca regulation and XB Cycle XB Cycle Combining computational models

23. Ca Ca Ca Ca Ca Ca Ca Must Use Monte Carlo Methods Relaxed Activated – Generates Force

24. Sample Results - Variation of Gamma Ca Ca n = 0 g Exponent on n = 1 n = 2

25. Some Results - F-Ca Relations Unpublished data from Dobesh et al., 2001, rat skinned fiber, SL = 2.15 mm

26. Mean Field Spatially Explicit Force Force Force Force 30x Residual force = incomplete relaxation [Ca] [Ca] 10x 10x Why correct F-Ca important? • [Ca] changes by ~10x while force changes by ~1000x in real muscle • Incorrect Force-Ca relationship leads to muscle that can not fully relax • In whole heart, muscle must relax to fill completely and pump efficiently 1000x

27. Need explicit consideration of spatial interactions Mass-action, mean-field model Spatially-explicit model Y. Yaniv, R. Sivan & A. Landesberg, Am J Physiol Heart Circ Physiol. 288, No. 1, H389-99. (2005). Rice et al., Pacific Symposium on Biocomputing, 2008

28. Mean Field Spatially Explicit Force Force Force Force 30x Residual force = incomplete relaxation [Ca] [Ca] 10x 10x Why correct F-Ca important? • [Ca] changes by ~10x while force changes by ~1000x in real muscle • Incorrect Force-Ca relationship leads to muscle that can not fully relax • In whole heart, muscle must relax to fill completely and pump efficiently 1000x

29. Combine 32 full sarcomeres into 1 myofibril which is a cell-level structure Moving to a terascale model of a myofibril 1 full 3-D sarcomere with hexagonal lattice Thick and thin filaments at molecular level

30. Decomposition for parallel processing

31. Decomposition scheme

32. Processing steps for each Dt If no D XB attachment Build A and K matrices in distributed manner with PETSc Update Markov states in each partition If D XB attachment Update locations by solving A.X=K with PETSc

33. Relaxation is not homogeneous Relaxation of myofibril showing inhomogeneity and spontaneous oscillations (unpublished data from B. Iorga and R. Stehle, 2006). time R. Stehle, M. Krüger, and G. Pfitzer Biophys J, October 2002, p. 2152-2161, Vol. 83, No. 4

34. Future directions - Whole heart model with large-scale parallelization on Blue Gene(collaboration with Universität Karlsruhe)

35. Spatial resolution is still too low ~ 200 elements 0.1 mm per cell 1010 elements

36. Blue Gene Machine - IBM T. J. Watson Research Center

37. Full heart and ECG calculation is not feasible F. H. Netter. Thieme. 1990 0.1 mm per cell 1010 elements Keller et al., 2007

38. Cell Models High degree of complexity Computation of potentials and concetrations Long simulation times Cellular Automaton Rule–based Pre–calculated potentials Short simulations times Modeling Electrophysiology Courtemanche et al. Cardiovasc Res. 1998

39. Experimental ventricular wedge preparation Yan, Shimizu, and Antzelevitch, 1998 Shimizu and Antzelevitch, 1998

40. Constructing a wedge model: Electrophysiology Clancy and Rudy, 2005 Silva and Rudy, 2005

41. Constructing a wedge model: 1-D cable and computing the ECG 100 mm - + 192 191 190 1 3 2 Endocardium Epicardium Transumural ECG

42. Sample ECG signals for the WT case M Cell Endo Epi

43. WT vs. Mut A vs. Mut B versions of IKs

44. Mut A Transmural ECG of Mut A and Mut B M Cell M Cell Endo Endo Epi Epi QT interval: 224ms(WT) 236ms(Mut A)

45. Next step: Move to 3-D wedge model Epi M-Cell Endo Torsade de pointe • 50x50x100 elements (=250 K) • Slow conduction velocity • Realistic size = 32 M elements Carruth and Silverman, 1980

46. Current approaches do not scale to large numbers of processors Kharche et al. (2008) in press Plank et al. (2006)

47. Decomposition for Blue Gene computer 1 0 3 2 1 0 3 2

48. Communication Framework Send/ReceivePointer Common Faces

49. Initialization Compute Vm Cycle 1 - Vm monodomain Compute I(Vm) bidomain Compute e Cycle 2 - e Compute I(e) CommunicationCycle

50. Results show good scaling to 8K CPUs