MUSICO Computational Platform
Welcome to the MUSICO PLATFORM for multiscale modeling of structure and function of the heart. This is an international effort of scientists developing computational models and quantitative, analytical tools for muscle research, closely tied with experiments designed to produce data that can be used to test model validity, flexibility, sensitivity, robustness.
This work will facilitate advancement of medical research in heart and skeletal muscle diseases, become a useful clinical tool for disease prediction and testing of novel therapies, and strengthen the fabric of international collaboration in the research community.
General multiscale model of muscle contraction includes two inter-connected models:
-
MUSICO Platform bridging the scales from inter-atomic (~ A) whole cell (~100-500 µm)
-
MUSICO Multiscale FEA and DSI bridging the meso-scales from molecular interactions (~ nm) to tissue (~ mm) or organ scale (tens of cm)
The multiscale problem of the heart. Structural data is defined from the scale of molecular bonds (10-10 m) and protein crystal structures to the 10 cm scale of the human heart. Time scales run from the pico-second molecular dynamics through the msec scale of protein-protein interactions to the 1 hz of the human heart beat. MUSCIO will integrate data from each level of organisation into a single comprehensive multiscale model.
Computational platform MUSICO (MUscle SImulation COde) provides realistic modelling of sarcomeric system with the aim to simulate a wide variety of experimental muscle behavior. The platform offers a modular program structure that allows extension and replacement of any part of sarcomeric system (calcium activation, cross-bridge cycle, sarcomere geometry, etc.). The current version of the MUSICO involves a number of sarcomere geometry models including the three-dimensional spatial models of multi-sarcomere geometry, multiple actomyosin cycle and calcium regulatory models. Nonlinear mechanical behavior of extensible filaments and crossbridges is addressed using iterative finite element scheme. Under development is inclusion of structural ancillary proteins MyBP-C, titin and nebulin/nebulette in 3D sarcomeric model.
The modular structure of the platform readily accommodates extension and replacement of any part of the sarcomere system ( activation, crossbridge cycles, sarcomere geometry, etc.) to incorporate new knowledge as it appears, allowing the platform to simulate muscle type or species specific differences. For example there are models of increasing complexity for the actomyosin cycle and regulation of contraction to account for increasing amounts (variables) of biochemical and mechanical (linear and nonlinear) information. Specific examples of complexity incorporated into MUSICO include the nonlinear mechanical behavior of extensible filaments and crossbridges, using an iterative finite element scheme, and the continuous flexible tropomyosin (Tm) chain. All of these systems are adaptable (for different muscle systems), scalable and parallelizable to address larger and larger ensembles as computing power permits. The current form includes the concept of distributed calculations of multi-scale domains in a mixed CPU–GPU environment.
The platform consists of multiple submodels:
Molecular Dynamics Simulations provide elastic characteristics of sub-molecular structural elements: in myosin S1, lever arm and neck region collectively composing crossbridge stiffness; in regulatory proteins tropomyosin, troponin and its subdomains, interaction of regulatory units with actin.
Molecular Kinetics and ATPase. The models of the experiments in solution: actomyosin ATPase cycle; - regulation of the ATPase cycle including & myosin induced changes in the thin filament structure, cooperativity of the filament, switching on and off of the activity. The models of thin filament regulation include Continuous Flexible Chain (CFC) and MacKillop-Geeves 3state cooperative model. The models are adaptable to include structural and kinetics changes between contractile protein isoforms or mutated proteins (myosins, regulatory proteins and ancillary proteins).
Models for Motility Assay Simulations. These models (under development) include regulation in motility assays environment incorporating geometrical factors associated with randomly oriented myosins on cover slips interacting with sliding actin filaments & simulations of different mutations, and mixtures of isoforms, and the effect cMyBP-C binding to actin or myosin head.
Multisarcomere 3D Explicit Stochastic Model. Multiscale modeling integrates the action of molecular (crossbridge) forces into a larger scale 3D sarcomere structure. A unified description of muscle function has been achieved by using a finite element model in which all geometrical factors and the constitutive elastic properties of filaments and crossbridges are all compounded in the system’s stiffness matrix. On the molecular scale, the kinetics of the underlying biochemical reactions defines the state rate transition matrix that integrates the biochemical reactions of the actomyosin cycle and its regulation, and the regulation of cell activity by undefined messengers or excitation-contraction-coupling pathways. Integrating these two aspects of the model, which operate on different length scales, and the activation and boundary conditions ( [ ] transients and external forces or displacements) requires iterative numerical procedures to predict the effects of strain-dependent modulation of the key biochemical reactions over multiple scales and nonlinear finite element analysis.
Sarcomere Structure and Explicit 3D Geometry. MUSICO simulations are defined by applying knowledge of the crossbridge (actin-myosin-ATP) cycle and its regulation in 3D lattice of striated muscle. The myosin (thick filament) and actin (thin) filaments are treated as elastically extensible, explicate positions of myosin heads (crossbridges) and associated myosin binding sites on actin filaments, and azimuthal position of regulatory chains.
Actomyosin Cycles. Actin-myosin interactions are defined in the context of 3D sarcomere geometry with discretely-positioned heads on myosin filaments and target zones on actin filaments. The strain-dependent actin-myosin kinetics are derived from reaction-energy landscapes So far this complex strain dependence is incorporated in 3 state cycle, 4 state cycle and 9 state cycle.
Thin Filament Regulation 3D Sarcomere Lattice. In current MUSICO platform continuous flexible chain model and MacKillop-Geeves 3state (cooperative) model are incorporated.
Thick Filament Regulation (planned). Thick filament regulation related to the super relaxed state (SRS) including stretch imposed activation is planned to be incorporated in the near future.
Structural Ancillary Proteins. In order to evaluate modulation of muscle contractility by ancillary proteins the modules incorporating titin, nebulin/nebulette and MyBP-C are currently under development.
X-ray Diffraction Patterns and Mechanical Data. Small-angle X-ray diffraction of contracting muscle during various mechanical protocols provides a view of sarcomere structure averaged over the area illuminated by the X-ray beam. Using MUSICO we assessed, for the first time, non-uniform changes in actin spacing, providing an estimate of the distribution of forces acting inside living cells. This provides a unique tool to measure the force acting on a single myosin or actin filament and therefore the possibility for more precise assessment of crossbridge forces.
Current stage of development of the above modules is shown in the Table.
MUSICO meso-scale models coupling FEA and DSI:
Simulating a heartbeat in whole heart is a complex, multi-scale problem. It is necessary to prescribe both, electrical and mechanical characteristics of the myocardium and it requires three principal components:
-
A multiscale model coupling molecular kinetics and its regulations with local material characters for macroscale continuous mechanics (finite element analysis).
-
Mesoscale structural organization of the heart tissue
-
Excitation-contraction coupling that is an essential step of electro-mechanical cardiac simulations
Coupling MUSICO with FE solver.
In nonlinear strain (displacement) based FE solvers, stress tensor in each integration point is calculated as a function of strains in the point, according to material model of the element. However, during muscle contraction, the material characteristic of muscle depends on local rate and history of deformation and cannot be measured and prescribed in integration points. Thus using multiscale formulation the instantaneous muscle characteristics should be calculated from the material model prescribed by a few kinetics and structural parameters. This can be achieved by coupling MUSICO and FE solver. For prescribed rate of deformation, i.e. shortening velocity of muscle fiber, MUSICO can calculate in any instant the local active stress (in direction of muscle fiber) and fibre stiffness. For known instantaneous muscle active stress and stiffness, along with passive material characteristics and current boundary conditions, FE solvers will provide updated muscle deformation and rate of deformation.
Multi-scale model of heat contraction. a) The 3D image of the hart discretized into FEs; b) The DTI tractography provides direction of the muscle fibres contained, shown for simplicity within a characteristic two-dimensional (2D) FE, including denoted integration points and the principal direction muscle fibers; c) elongation of an individual muscle fiber, at the indicated spatial scale and under stress; d) Cross-bridge kinetics mode.
To achieve this, at muscle tissue level detailed information of the orientation of the fibers is necessary to prescribe complex constitutive equations, as orthotropic models. We have developed methodology for determining orientation of the cardiac fibers at the organ level using diffusion spectrum imaging (DSI).
Determining Fiber Orientation for Cardiac Computational Modeling.
The essential information necessary for finite element analysis is the orientation of fibers due to two reasons: (i) the direction of active forces is aligned with the direction of fibres and (ii) the electricity travels much faster along muscular fibers than across. Thus, a correct anatomical description of the arrangement of muscle fibers is probably the most important element in obtaining a realistic model of the heart. Cardiac fibers are spatially organized to enable coordination to allow the contraction of heart cavities in a coordinated, stable and efficient way. Today we can achieve much better description of the fibers using a technique called diffusion tensor imaging, using magnetic resonance map of diffusion along biological fibers. The variations of this technique are diffusion spectrum imaging (DSI) and generalized Q-space imaging (GQI). We use DTI, DSI and GQI, depending on the sequence of acquired images, to validate the theory with experiments at organ level. We show in figure below the methodology how we can extract the direction of fibers from fiber tracks.
Determining fibre orientation across the ventrical walls. (A) Reconstructed fibres from DTI imaging of the heart; (B) at each cross section the fibres create sheets that change orientation across the thickness of the wall; (C) illustration of the distribution of angles across left ventrical wall (a) and the muscle fibres form layers 3 -4 cells thick (b). The layers are interconnected by collagen; (D) newly developed method to obtain direction of fibres at any position in arbitrary cross-section (yellow arrows). This directions are directly implemented in finite elements integration points along with other active and passive material characteristics; (E) reconstructed directions of muscle fibres across and along ventrical wall; (F) The reconstructed DTI images in concert with fine finite element mesh realistically represent the geometry of fibrous and layered structure of ventrical walls.
Architectural basis for impaired cardiac mechanics.
Considering the mechanism of the architectural phenotype associated with cMyBP-C3 mutations, we employed magnetic resonance imaging (MRI) methods to track molecular diffusion and diffusion sensitivities with tractography. Employing these methods, we found that normal hearts displayed a gradual variation of helix angle from sub-epicardium to sub-endocardium, similar to the helical fiber patterns portrayed by Streeter et al. Moreover, the ablation of cMyBP-C-/- is associated with loss of the predicted variation of transmural helix angle configuration and a loss of intervoxel orientational coherence in the region spanning the mid-wall to the endocardium. Depth-specific differences of myoarchitecture may impact on cardiac tissue mechanics in that cMyBP-C-/- hearts possess a shortened contraction time and achieve greater peak elastance at the onset of ejection compared to the wild type. Altered transition of helix angles may also enhance the extent of initial elastance by means of faster force activation. The more rapid onset of elastance and lack of helical fiber development from the mid-wall to the sub-endocardium promote a shorter contraction time. We show that cMyBP-C-/- mouse hearts display a distinct loss of fiber helicity, implying that certain mutations in cMyBP-C-/- are critical for the architectural “design” of the myocardium and may represent an important branch-point in the modulations of cardiac contractility.
Imaging the myoarchitecture of the ventricular wall with GQI
Architectural phenotype associated with constitutively impaired M domain phosphorylation
Modeling Cardiac Electro-mechanical Coupling.
The cardiac computational model requires solution of three components: the electrical signal, the mechanical deformation and the excitation-contraction coupling that links both problems. The equations of the complex models describing the heart can be only solved on a computer using numerical techniques. The methods are based on a finite element approach (FEM), focusing on exploiting a high parallelization software to be solved in high performance computers with high efficiency. The time evolution of the system is calculated with finite differences (FD), using an Euler scheme.
The excitation-contraction coupling is the phenomenon by which the fibers contract after the electrical activation propagates through the myocardium. Thus, the contraction is triggered by electrical depolarization, which is initiated at specialized myocyte cells, called Purkinje cells. The Purkinje network is part of the cardiac conduction system that is responsible for the fast and coordinated distribution of electrical impulses in the ventricles. Purkinje system is electrically connected to the myocardium muscle through transitional cells or Purkinje-myocardium junctions (PMJ) distributed throughout the ventricular sub-endocardial layer. The myocardium is generally regarded as an anisotropic material with respect to its fiber orientations, and there is strong evidence of its orthotropic behavior. Cardiac fibers’ orientation requires a vector field obtained from either histological studies made in ex-vivo hearts, benchmark synthetic models or diffusion spectrum (magnetic resonance) imaging (DSI).
The use of DSI has become the gold standard measure for fiber orientation. However it is not yet generally available for in-vivo human studies due to the lack of resolution and quality. Sometimes, warping of fiber orientations from models of alternate species have been used to fill in the gap of information. Only recently, there has been a thorough study where significant differences in the myocardial fiber structure as represented by the fiber helix angle were observed between species like the mouse, rabbit, sheep and humans. The effect of fiber orientation on the propagation of the electric wave has been assessed at organ level by reaction-diffusion systems. This system reproduces continuous approximations of the excitation of cardiac muscle. In collaborative effort with Alya System, i.e. Alya Red project of the Barcelona Supercomputing Center for simulating a human heart, we plan to develop complete model of the heart including the above multiscale FEM, DSI, and electrophysiological activation.