Feedback, lineages and cancer; John Lowengrub, University of California at Irvine, USA

Most tissues are hierarchically organized into lineages. Tumors arise when the carefully regulated balance of cell proliferation and programmed cell death (apoptosis) that ordinarily exists in normal homeostatic tissues breaks down. Traditionally, cancer cells are assumed to acquire a common set of traits. However, not all proliferating cells in a tumor matter seem to matter equally and tumor cells apparently progress through lineage stages regulated by feedback. The reasons for this are still not well established. In this talk, we will present new mathematical models to simulate the spatiotemporal dynamics of cell lineages in solid tumors. We account for protein signaling factors produced by the cells and nutrients in the microenvironment. Our models demonstrate the development of heterogenous cell distributions and the formation of niche-like environments for stem cells. We demonstrate that feedback processes play a critical role in tumor initiation, progression, the development of morphological instability, and the response to treatment.

HPC & The Inventor’s Dilemma: Abandoing 20 years of Molecular Simulation Knowledge and Moving to GPUs; Eric Lindahl, SciLifeLab Stockholm, KTH, Sweden

Molecular dynamics simulations are undergoing a second revolution. Computers are suddenly fast, large, and cheap enough that it is possible to perform simulate a large number of biomolecular processes that were completely out of reach only a few years ago, for instance the opening or closing of an ion channel. However, the free lunch in terms of constantly faster processors is likely over – while computers get cheaper per processor, the improved performance now has to come from increasing parallelism, which is a huge problem for a computational algorithm such as molecular dynamics that historically has not scaled well for small problem sizes. Second, after 20 years of work in a field programs tend to be large and complex, which makes them difficult to turn around. Here, I will describe our early efforts on high-end optimization & parallelization of a scientific code, and in particular how we decided to revisit the entire problem from scratch a few years ago. Not only did this lead to a molecular simulation implementation that is by far the fastest of any on GPUs, but surprisingly we could take many of the ideas back to modern CPUs and improve scaling considerably.

HPC-ABGEM: High-performance agent-based general ecosystems modeling; Simone Callegari, John D. Weissmann, Natalie Tkachenko, Wesley P. Petersen, George Lake, Christoph P. E. Zollikofer (ETH and University of Zurich)

The HPC-ABGEM project within PASC aims to develop an agent-based simulation code and model library suited for HPC. Agent-based modeling is a powerful numerical approach, in which emergent large-scale properties of a complex system are simulated bottom-up, starting from the behavioral rules of its building blocks (the "agents"). This approach is used in many fields including ecology, sociology, and economics. It is well suited in particular to model population and ecosystem dynamics across scales, from individuals to species, from fractions of a lifetime to evolutionary timescales. Here we present some of the computational challenges of agent-based modeling at very large scales, describe our HPC approach, and show preliminary results from tests of stochastic population dynamics. Once ready, this code will be used to study the world-wide out-of-Africa dispersal of anatomically modern humans, and the interaction between Homo sapiens and Neanderthals in Europe during the Late Pleistocene.

Bone structure analysis on multiple GPGPUs; Peter Arbenz, Cyril Flaig, Daniel Kellenberger (ETH Zurich)

Osteoporosis is a disease that affects a growing number of people by increasing the fragility of their bones. The lifetime risk for osteoporotic fractures in women is estimated close to 40%, for men risk is 13%. Since global parameters like bone density do not admit to predict the fracture risk, patients have to be treated in a more individual way. Today’s approach consists of combining 3D high-resolution CT scans of individual bones with a micro-finite element analysis. To improve the understanding of bone, large scale computer simulations are applied.

In recent years we have developed a fast, scalable and memory efficient solver for such problems called ParOSol. The equations of linear elasticity in pure displacement formulation are discretized by the finite element method. The resulting linear system of equations is solved by the conjugate gradient algorithm with a multigrid preconditioner. The solver is matrix free on all grid levels. Note that the principle obstacle in this problem is the very complicated computational domain entailed by osteoporotic bones.

We have ported the solver to GPGPUs to profit from their exorbitant compute capabilities. We discuss how we modified the CPU code to run successfully on multi-GPGPUs. We present the implementation of ParOSol on Tödi. The code shows perfect weak scaling up to 256 GPGPUs. The largest problem solved on the machine contained about 8 billion voxel elements corresponding to about 25 billion degrees of freedom. Compared to the pure MPI code run on the 16-core nodes of Tödi we observe speedups of the MPI-Cuda code beyond five.

Computational simulation of heart function with an orthotropic active strain model of electromechanics; Toni Lassila, Simone Rossi, Alfio Quarteroni (EPFL), Ricarco Ruiz-Baier (UNIL)

Fully integrated models of total heart function including the fluid dynamics inside the ventricle are on the horizon as computational algorithms become more efficient, allowing more comprehensive simulations of cardiovascular pathologies. Due to the multiscale and multiphysics nature of heart activity, a careful consideration between different algorithmic choices (monolithic vs. segregated, tightly vs. loosely coupled, explicit vs. implicit) has to be made in order to obtain the best compromise between stability and scalability. Our aim is to construct a baseline multiphysics integrated heart model that can be used to study the interaction between four basic fields related to the heart function: the cell-level electrophysiology, the subcellular activation mechanisms, the elastic deformation of the tissue, and (optionally) the ventricular fluid mechanics.

To model the electromechanical muscle contraction we use a stretch-dependent activation model based on thermodynamically consistent transversely anisotropic active strain model that is able to capture correctly the transmural thickening of the tissue. The passive mechanics model is adapted from the seminal work of Holzapfel and Ogden and includes both fiber and sheet directional orthotropy. For the ventricular fluid-solid coupling we opt for a monolithic approach as it provides in our experience the best stability and computational efficiency for large displacement fluid-structure interaction problems. Additional ongoing work is related to the modelling of heart valves and their effect on ventricular fluid dynamics, in order to study e.g. the effect of valve pathologies on the mechanical work performed by the ventricles.

A computational model for total heart function has been implemented in the open-source finite element library LifeV (http://www.lifev.org) and recently released to the public. It is based on parallel solvers for the electrophysiology, solid mechanics and fluid mechanics subproblems. Mesh partitioning is performed with ParMETIS, and the linear problems are solution strategy based on multigrid or algebraic additive Schwarz preconditioners implemented within the Trilinos library.

Red blood cells and related functions of the human spleen; Igor V. Pivkin (Universita della Svizzera italiana), Zhangli Peng (MIT), Ming Dao (MIT)

The human spleen plays a key role in innate and adaptive immune responses, in selective clearance of aged or abnormal red blood cells (RBCs), and in removal of pathogens present in the blood. It is also central to the pathophysiology of several frequent and potentially severe diseases such as inherited red cell membrane disorders, hemolytic anemias, and malaria. We will present the results of simulations of RBC motion through the smallest openings in the human spleen. The modeling approach is based on the Dissipative Particle Dynamics (DPD) method.

LBIBCell: A Cell-Based Simulation Environment for Morphogenetic Problems ; Simon Tanaka, Prof. Dagmar Iber (D-BSSE, ETH Zurich)

We developed a cell-based simulation framework to study morphogenetic problems. The tissue model is based on the IBCell model approach, which models interactions between the viscous interstitial fluid and cytoplasm, and the elastic structures such as cell membrane, cell-cell junctions and cytoskeleton using the immersed boundary method. The tissue model is coupled to a system of advection-diffusion-reaction partial differential equations to model signaling, and thus is predestinated to study cell-based mechanisms based on the interaction of tissue dynamics and signaling, as appearing in many morphogenetic problems.

We studied the dynamics of a Turing-type reaction-diffusion signaling system, which locally controls the proliferation rate, and demonstrate that, in comparison to the continuous formulation conterpart, surprising effects (such as dividing cell clusters) emerge.

The software will be published under the GPL license.

SEM++: a particle model of cellular mechanics, growth, signaling and migration; Gerardo Tauriello (ETH Zurich), Florian Milde (ETH Zurich), Petros Koumoutsakos (ETH Zurich)

We present a discrete particle method to model biological processes from the sub-cellular to the inter-cellular level.

Compared to continuum models of biological tissues, a discrete model allows the exploration of single cell dynamics that are important for cell division patterns, growth and adhesion. If we model cells at a continuum level, a single cell does not exist. Therefore, the emergence of phenomena initiated by single cells cannot be studied in those models.

Our discrete model is an extension of the subcellular element method which represents each cell as a collection of particles. Each particle represents a coarse grained portion of the cytoskeleton which is embedded in a viscous fluid—the cytoplasm. Particles interact through a parametrized force field to model cell mechanical properties. We extend the model by including distinct particles for the nucleus and the cell membrane. We include hydrodynamic interactions to capture the fluid mediated forces within the cytoplasm. We introduce a model for cell growth and proliferation driven by the cell-cycle and we account for actin polymerization based migration. We propose extensions to model cell-cell mediated signaling processes that can determine the cell phenotype.

We discuss the guiding design principles for the selection of the force field and the validation of the particle model using experimental data. The proposed method is integrated into a multiscale particle framework for the simulation of biological systems.

Exploring the unfolded ensemble of proteins; Albert Ardevol (ETH Zürich), Michele Parrinello

Protein flexibility plays an essential role in life. In the last few years, it has become apparent that biologically functional processes cannot be understood unless the dynamical behavior of the system is accounted for. For instance, protein flexibility has been found of outmost importance in several processes such as ligand binding, signal transduction and catalysis, to name but a few. This ensued a change in the structural biology paradigm. Theories like allosteric binding or induced-fit theory were proven, and other like conformational selection, conformational shift mechanism or barrier-less downhill folding have arisen. From a theoretical perspective, the thermodynamic and kinetic behavior of flexible biomolecules have has been successfully addressed by the investigation of their underlying free energy landscape (FEL). However, these FELs are usually very complex, with multiple metastable states and transitions that exceed the microsecond time scale and even the space of collective variables that define unequivocally all relevant structures and motions is not straightforward. In the present work, we applied a combination of two methods recently developed in the group (namely the well-tempered ensemble[1] and the sketchmap[2]) to study the conformational ensemble of a short peptide that folds into a beta-hairpin conformation at room temperature, yet putting the focus on the unfolded state of the peptide. We show how the well-tempered ensemble and metadynamics can be used to enhance the performance of a parallel tempering simulation to sample complex conformational landscapes. The tertiary structures explored by the peptide during the simulation were analyzed using sketchmap to transform the multidimensional landscape defined by all the backbone torsional angles onto a two dimensional surface, in which the relative free energies were calculated. We show how, in addition to a random coil, the unfolded state is composed by several metastable basins with a defined structure that are in fast equilibria with the folded beta-hairpin. Finally, we analyzed the effect that certain mutations have on the relative free energies of the folded, misfolded and unfolded conformations.

[1] Bonomi, M.; Parrinello, M. PRL 2010, 104, 190691-4 [2] Ceriotti, M.; Tribello, G.; Parrinello, M. PNAS 2011, 108(32), 13023-28

Protein-ligand unbinding kinetics and pathways through metadynamics; Pratyush Tiwary (ETH Zurich and Università della Svizzera italiana), Vittorio Limongelli (University of Naples), Matteo Salvalalgio (ETH Zurich and Università della Svizzera italiana), Michele Parrinello (ETH Zurich and Università della Svizzera italiana)

Metadynamics is a commonly used and successful enhanced sampling method in atomistic simulations, with applications that range from nanotechnology to biophysics. By the introduction of a history dependent bias which depends on a restricted number of collective variables it can explore complex free energy surfaces characterized by several metastable states separated by large free energy barriers. Recently we extended its scope by introducing a simple yet powerful method for calculating the rates of transition between different metastable states. A statistical analysis of the distribution of transition times allows us to verify a posteriori the accuracy of our calculations. We apply our method to predict the average escape time for a ligand unbinding from a protein, which is a very important parameter in a drug design strategy. Specifically, we study the the trypsin-benzamidine system in explicit solvent and calculate the average escape time, lifetimes of stable intermediates and fluxes for various unbinding pathways.

Improving I/O Performance of ClustalW-MPI Multiple Sequence Alignment Software; Soon-Heum Ko (Linkoping University), Plamenka Borovska, Veska Gancheva (University of Sofia)

We introduce the deployment of a parallel I/O technique to an MPI version of ClustalW multiple sequence alignment software. This work aims to accelerate the I/O performance of ClustalW-MPI code so that the I/O part of this code can fit to aligning massive numbers of biological sequence chunks whose data size can range to a Gigabyte level in total. The work has started by understanding the structure of the code and I/O routines. This code adopts a master-slave parallelism so that the master performs all global comparisons and dynamically distributes the embarrassingly parallel task to each slave. The entire sequence input is first stored in a master and passed to each slave by using the point-to-point communication operation (MPI_Send), which is expected to benefit by deploying MPI-I/O functions. On the other hand, the output operation is hard to be parallelized since it is preceded by a global comparison of alignment results at a master processor. Thus, the main effort has been put towards employing the parallel I/O technique to reading the sequence input and tuning this implementation. The parallel I/O function has been implemented in a way that a subset of MPI cores simultaneously performs the I/O operation and locally broadcast to individual neighbours so that the code is less sensitive to the stability of the parallel file system. We performed an intensive benchmark to tune important parameters of (1) number of I/O ranks, (2) read data size per each MPI_Read instruction, and (3) tunable parameters provided by MPI-I/O library. In our experiments whose file sizes range from tens of kilobytes to one gigabytes and MPI sizes range from a thousand to four thousand cores, the optimal number of I/O ranks was either 4 or 8. Letting all MPI ranks to simultaneously conduct the read operation is highly avoided in that it reduces the file read speed and increases the risk of encountering out-of-memory error in the MPI-I/O library. The read array size is recommended to be as big as the file size so as to minimize calling MPI_Read function. On the other hand, increasing this read array size also increase the risk of hitting out-of-memory error from MPI-I/O library. With these conditions applied, parallel I/O outperforms the serial I/O by 2 – 4 times at tens or hundreds of megabyte datasets. Out of tunable MPI-I/O parameters, we examined the performance variation by changing cb_nodes (which denotes the maximum number of aggregators) and cb_buffer_sizes (which expresses the internal buffer size to transfer at a single instruction). The effect of cb_nodes was not observed because it is much dependent of hardware structure. Regarding cb_buffer_sizes, a smaller size reported the better performance, which is different from our common understanding that a larger buffer size will outperform because it reduces the number of read operation. It is hard to argue that this is the general characteristics, but we can affirm that the small cb_buffer_size is preferred for some applications which inputs the character strings. Finally, the current MPI-I/O implementation has been employed on ClustalW-MPI code. Experimented sequence data size is 16 megabyte and it contains 8,999 sequences. Experiments have been performed on a BlueGene/Q system and the number of cores varies from 512 to 8192. The performance with parallel I/O outperforms the original implementation in all experiments and the difference widens as the number of cores increase so that parallel I/O finally approaches 6.8 times faster than the original operation. It concludes us that the current parallel I/O interface design is reasonable for application to many bio-informatics codes of similar kinds.

OmicABEL: Story of a successful interdisciplinary collaboration; Diego Fabregat-Traver (AICES, RWTH Aachen, Germany), Yurii S. Aulchenko (SD RAS, Novosibirsk, Russia), Paolo Bientinesi (AICES, RWTH Aachen, Germany)

We discuss the main hurdles and results of a four-year long collaboration between a group of computer scientists and a group of computational biologists. Despite the difficulties inherent to every interdisciplinary collaboration, we can safely affirm that ours is a story of success: we recently released OmicABEL, a high-performance library for large-scale genome-wide association studies (GWAS) that attains remarkable speedups over the state-of-the-art tools, and most importantly, we attracted more research groups and paved the way for further joint research.

Due to the boost in the amount of available genomic data, computational biologists are eager to perform GWAS analyses of ever increasing size: a typical large analysis involves the computation of billions of generalized least-squares (GLS) problems and the processing of terabytes of data. Unfortunately, neither general nor domain-specific libraries provide a viable solution, requiring unreasonable computing resources and time to completion. Our main goal thus has been the development of efficient and scalable algorithms and routines for the computation of large-scale GWAS.

The main limitation in available tools is that they are based on the traditional black-box approach where the individual GLSs (or small subsets of them) are addressed in isolation. Instead, the key to a viable solver is to consider the entire grid of problems as a whole, taking advantage of the specific correlation among them, and to uncover and exploit as much domain knowledge as possible (e.g., the special structure of matrices). Based on this principle, our application-tailored algorithms lower the computational cost by orders of magnitude. Combining the efficiency of our algorithms with high-performance techniques to attain scalability, and out-of-core techniques to effectively manage datasets residing on disk, our library outperforms the state-of-the-art software by factors of up to 1000. Analyses that previously required the use of supercomputers can now be performed on one multi-core node in a matter of hours.

While a remarkable result, the outcome of our collaboration goes beyond OmicABEL. The collaboration turned into a breeding ground, with multiple research groups involved in the development of a variety of projects. Finally, we highlight the most important lesson we learned from this collaboration: There is a huge disconnection between the field of high-performance computing and the application domains. When both worlds come together, however, the quality of the outcome grows exponentially. In order to advance in the field of simulation sciences, it is imperative that experts from different fields sit together and work on filling this gap. Close collaboration and constant communication are crucial to overcome language barriers and to achieve satisfactory results.