Description Language of Calculation Scheme for Automatic Simulation Code Generation

1. Description Language of Calculation Scheme for Automatic Simulation Code Generation Akira Amano College of Lifescience, Ritsumeikan University CellML workshop 2012

2. HD Physiology Project • Project Leader: Prof. Y. Kurachi (Osaka Univ.) • 2010-2015 • http://hd-physiology.jp/ CellML workshop 2012

3. HD Physiology Project • Japanese title of this project “Integrative Multi-level Systems Biology for …” • target • Multiscale in Biology • construct comprehensive simulation model of drug induced heart arrhythmia CellML workshop 2012

4. HD Physiology Project • groups • platform (Kitano (SBML), Asai (isML PHML)) • cardiac biology • pharmacokinetics • integrate models of ... • whole heart excitation propagation • cell model including drug effect • circulation model • pharmacokinetics / pharmacodynamics CellML workshop 2012

5. Characteristics of Biological Function Models (from IT view) • Biologicalfunction simulations are different from conventional models in: • The complexity of the models • Multiple scale and physics • (nano-sec, hour, Chemistry, dynamics, fluid ...) CellML workshop 2012

6. Characteristics of Biological Function Models (from IT view) • The complexity of the model • Multiple scale and physics cardiac cell model latest model number of equations： more than 400 number of variables： more than 400 ion channels ion transporters contraction mechanism energy metabolism CellML workshop 2012

7. Characteristics of Biological Function Models • The complexity of the model • Multiple scale and physics Simulation of the heart beat fluid mechanics structural mechanics Coupling Calculation electrophysiology chemical reaction Coupling Calculation Scheme is usually built in in the simulation tools. CellML workshop 2012

8. Software Crisis! in 1970s: Software size increased exponentially. Large scale software project failed. Merlin Dorfman ”Software Engineering” IEEE,1997. Software Engineering Structured Programming Languages Integrated Software Tools CellML workshop 2012

9. Software Engineering Tools for Complex Biological Models Efficient handling of complex models ⇒ use of Markup Languages Cell physiological model description Field (Organ) model description CellML FieldML Systems Biology model description PHML SBML Declarative Description Models CellML workshop 2012

10. Software Engineering Tools for Complex Biological Models Declarative description: every equation is true at any time not declarative declarative • High flexibility • from mathematical description • combination of models, boundary cond. • Avoiding bugs • by automatic processing of models CellML workshop 2012

11. What to be supported • Lumped parameter models (ODE) simulation tools and markup languages CellML, ... calculation methods are simple (usually builtin) • Distributed parameter models (PDE) simulation tools and markup languages (complex) FieldML?, ... calculation methods are complex no general simulation software? • Coupling simulation models simulation software are developed manually Our approach: propose description language for calculation methods CellML workshop 2012

12. What to be supported • Method: • Designing declarative description of • coupling calculation schemes. • use recurrence relations • start from single model discretization • next consider coupling calculation CellML workshop 2012

13. Single Model Input Equations • Lumped parameter models (ODE) • discretize in time • Euler, Runge Kutta, ... • Distributed parameter models (PDE) • discretize in space (and time) • center difference, upwind difference, ... CellML workshop 2012

14. Discretization FTCS (forward time center space) original model equation discretized equation CellML workshop 2012

15. Discretization Crank-Nicolson original model equation discretized equation • calculation method becomes complex • we can separate discretization scheme from calculation method CellML workshop 2012

16. Process of combining Model Equation with Calculation Scheme CellML Compiler available at sourceforge CellML/ PHML Model Equations Recurrence Relation Equations Executable code equation to program easy! CPU, GPU PC cluster, super computer (K?) low efficiency but high and safe productivity TecML Calculation Scheme Relation (RelML) 1st step 2ndstep Variable Correspondence information between two files CellML workshop 2012

17. Declarative Description of Calculation Schemes • Single Model Calculation schemes: • mathematical relation between discretized variables • calculation method of unknown variables • TecML • (Time Evolution Calculation Markup Language): • Description language to describe • time evolution calculation method • Euler, RungeKutta … • spatial difference scheme • center difference, upwind difference, ... CellML workshop 2012

18. Concept of Calculation Scheme Description by Recurrence Relation relation model equations Recurrencerelation RelML time discretization TecML (declarative description) 1st step: generate recurrence relation from model equations CellML workshop 2012

19. Concept of Calculation Scheme Description by Recurrence Relation Recurrencerelation Executable code x[0] = xinit; tn=0; do{ r[tn] = x[tn]*x[tn]*x[tn]; k[tn] = - r[tn]; x[tn+1] = x[tn] + k[tn] * delta; tn = tn + 1; } while(!((tn-1) == 1000)); xend = x[tn-1]; (declarative description) 2nd step: generate program to calculate recurrence relations CellML workshop 2012

20. Experiments on FHN, LR-I CellML workshop 2012

21. Coupling Calculation Schemes • Coupling calculation schemes: • also mathematical relation between discretized variables in several models • calculation method of unknown variables • TecML • Description language to describe • coupling calculation scheme • cell model + mechanical model • convergence calculation • different time step models CellML workshop 2012

22. Example: Double Time Loop Physiologically based Pharmacokinetic model calculation time order: hour time order: millisec CellML workshop 2012

23. Example: Double Time Loop do { dose concentration calculate whole body model with large time step do { concentration cell state calculate each cell model with small time step } while ( !time1end() ); } while ( !time2end() ); CellML workshop 2012

24. Example: Double Time Loop Recurrencerelation Executable code c[0] = 1.0; tn=0; do { s[tn][0] = s[tn-1]; tm = 0; do { ks[tn][tm] = f2(c[tn],s[tn][tm]); s[tn][tm+1] = s[tn][tm] + ks[tn][tm]*delta1; tm++; } while (!(tm-1)==1000); s[tn] = s[tn][tm]; kc[tn] = f1(c[tn],s[tn]); c[tn+1] = c[tn] + kc[tn]*delta2; tn = tn + 1; } while(!((tn-1) == 1000)); cend = c[tn-1]; inner while loop outer while loop CellML workshop 2012

25. Example: Convergence Calculation The dynamic model of the simple heart CellML1:Cardiac Cell Model CellML2:Blood Circulation Model CellML3:Heart Finite Element (FE) Model CellML workshop 2012

26. Example: Convergence Calculation CellML1:Cardiac Cell Model CellML1: Cardiac Cell Model dT/dt = kT = f1 (L, T) CellML2:Blood Circulation Model CellML2: Blood Circulation Model dV/dt = kV = f3 (P) CellML workshop 2012

27. Example: Convergence Calculation CellML3:Heart Finite Element (FE) Model CellML3: Heart FE Model FE model can only calculate volume from pressure. PLV is calculated by searching. V= f4 (P, T) L = f5 (V) nP = fPnext (P, V, Vtarget, P1, P2) nP1 = P nP2 = P1 e = |V – Vtarget| fpnext(P1,P2) P1 P2 Vnext Vtarget V1 CellML workshop 2012

28. Recurrence relation CellML workshop 2012

29. T[0]=Tinit; V[0]=Vinit; P[0]=Pinit; tn = 0; do{ Pr[tn][0]=P[tn]; P1r[tn][0]=P[tn]; P2r[tn][0]=P[tn]; tm = 0; do{ Vr[tn][tm]=f4(Pr[tn][tm],T[tn]); Lr[tn][tm]=f5(Vr[tn][tm]); nPr1[tn][tm]=fPnext((Pr[tn][tm],Vr[tn][tm], V[tn],P1r[tn][tm],P2r[tn][tm]); nP1r1[tn][tm]=Pr[tn][tm]; nP2r1[tn][tm]=P1r[tn][tm]; er[tn][tm]=|Vr[tn][tm]-V[tn]|; Pr[tn][tm+1]=nPr1; P1r[tn][tm+1]=nP1r1; P2r[tn][tm+1]=nP2r1; tm = tm + 1; }while(!((er[tn][tm-1])<epsilon)); P[tn]=Pr[tn][tm-1]; L[tn]=L[rtn][tm-1]; dTdt0[tn]=f1(L[tn],T[tn]); dVdt0[tn]=f3(P[tn]); T[tn+1]=T[tn]+(dTdt0[tn]*delta); V[tn+1]=V[tn]+(dVdt0[tn]*delta); tn = tn + 1; }while(!((tn-1)==1000)); Tend=T[tn-1]; Vend=V[tn-1]: Executable code CellML workshop 2012

30. 2nd step: calculation method • Coupling calculation schemes: • also mathematical relation between discretized variables of several models • calculation method of unknown variables • current status: • order of calculating eqns are analyzed • partly DAE calculation is supported • (numerical differential method) • simple GPU support CellML workshop 2012

31. advantages of our approach • use of declarative description is flexible in ... • combination of various models • combining with boundary conditions • two step approach • various target language / hardwares / libraries • language: C, C++, Java, Java BigDecimal • hardware: GPU, super computers • library: various math library CellML workshop 2012

32. Example: parallel code (GPGPU) Single CPU CellML Simulation Program Generator C language GPU CUDA C 2D excitation propagation (LR1) CPU (core i7 880) GPU (Tesla C2050) GPU was 60 times faster than single CPU We can use any cell model to generate GPGPU 2D excitation propagation program. The 1st HD Physiology International Symposium 20th Jan, 2012

33. current project webpage • available at sourceforge (http://sourceforge.net/projects/cellmlcompiler/) • only single model time discretization available • Euler, Runge Kutta (explicit method only) • fixed geometry excitation propagation with GPU CellML workshop 2012

34. Summary • Declarative description models are useful for simulation with complex models. • CellML, SBML, FieldML, PHML, ... • Declarative description of coupling calculation scheme • proposed: TecML • With declarative description coupling calculation scheme, simulation program can be generated with high flexibility. • various boundary conditions • parallel code for super computers CellML workshop 2012