Simulating and modelling complex biological systems in computational life sciences requires specialized software tools that can perform medical image data-based modelling, jointly visualize the data and computational results, and handle large, complex, realistic and often noisy anatomical models. The required novel solvers must provide the power to model the physics, biology and physiology of living tissue within the full complexity of the human anatomy (e.g. neuronal activity, perfusion and ultrasound propagation). A multi-physics simulation platform satisfying these requirements has been developed for applications including device development and optimization, safety assessment, basic research, and treatment planning. This simulation platform consists of detailed, parametrized anatomical models, a segmentation and meshing tool, a wide range of solvers and optimizers, a framework for the rapid development of specialized and parallelized finite element method solvers, a visualization toolkit-based visualization engine, a Python scripting interface for customized applications, a coupling framework, and more. Core components are cross-platform compatible and use open formats. Several examples of applications are presented: hyperthermia cancer treatment planning, tumour growth modelling, evaluating the magneto-haemodynamic effect as a biomarker and physics-based morphing of anatomical models.
Computational life sciences (CLS) have special requirements that are not sufficiently well satisfied by existing commercial simulation platforms. These include for example:
— the ability to perform simulations based on medical image data (to define the model geometry, boundary conditions and material parameter distributions);
— the ability to handle large, complex, realistic and often noisy anatomical models (e.g. generation of high-quality meshes based on segmentation label fields);
— the ability to jointly visualize and post-process simulation results and measurement data (e.g. imaging data);
— novel solvers providing the power to efficiently simulate the physics of living tissue (e.g. perfusion and thermoregulation models in tissue heating simulations, acoustic propagation in focused ultrasound modelling) within the full complexity of the human anatomy; and
— novel, specialized solvers for biological/physiological processes (e.g. EM–neuron interactions and tissue/tumour growth models)
The solvers have to be high-performance computing (HPC) enabled to be able to handle such complex models within an acceptable time.
To illustrate how a single device can require the modelling of various physical and physiological aspects as well as their interplay, consider a stent: the stent influences blood flow (fluid dynamics) and needs to be deployed (mechanical problem including contact). The flow results in forces and stress at the vessel wall (fluid–solid interaction), which in turn can be related to rupture risk or vessel wall remodelling (complex multi-physics–multi-scale problem involving for example coupled convection–reaction–diffusion (CRD) problems). Additionally, the impact of a stent on clotting or on vascular degenerative diseases, such as aneurysms, can also be of interest (multi-physics–multi-scale modelling coupling flow, solid mechanics and mass transfer). Safety considerations in the presence of high fields, e.g. during magnetic resonance imaging (MRI), include radiofrequency (RF)-induced heating (electro-dynamics and thermal), induced forces and torque, as well as nerve stimulation (neuron dynamics). Other safety aspects are related to, for example, mechanical wear and tribology. The choice of the stent type and its deployment strategy can require patient-specific, image-based treatment planning (segmentation, model generation). In the case of drug-eluting stents, there are additional aspects related to drug transport and activity (CRD). Modelling is useful when designing and optimizing stents, for clinical application and for investigating the underlying mechanisms of their physical and biological interaction with human physiology. A single device thus requires (multi-)physics/biology/physiology modelling on multiple scales, often involving anatomical models and (patient-specific) imaging-based information that rely on powerful solvers capable of handling complex problems.
In this paper, a multi-physics simulation platform developed and optimized for the needs of CLS and medical technology is presented. It satisfies the requirements detailed above and overcomes many of the the limitations of existing commercial software with regard to image-based modelling, simulations involving complex anatomical models, HPC-enabled solvers optimized for physical effects in living tissue, physiology/biology modelling, and joint visualization of simulation, imaging as well as measurement data. In this paper, its application is illustrated with a series of examples highlighting some of this functionality.
2. Material and methods
The simulation platform covers all components from medical image segmentation to the provision of detailed anatomical models, (parametrized) geometrical modelling, simulation set-up, domain discretization, solving (including management of hardware resources and coupling of solvers), optimization, post-processing, visualization, scripting and application development.
Fundamental to the framework concept is the rigorous modularity, the cross-platform compatibility of core components, and the strong integration of Python scripting  and the visualization toolkit (VTK) . In addition to graphical user interface (GUI)-based interactions, all functionalities are interfaced such that they can be accessed at a high level through Python scripts. Major components are encapsulated as VTK filters. This allows these components to be used as parts of VTK pipelines, and data to be passed in VTK formats (as well as in other open formats, such as the eXtensible Data Model and Format XDMF ), post-processed, and visualized with the full VTK functionality enhanced by a self-developed OpenGL visualization engine. A Python-interfaced GUI software development kit (SDK) allows the flexible extension of the provided GUI with user-defined plug-ins or the implementation of own applications by means of building blocks from the platform. Modelling can thus be performed within the GUI, through user-generated Python scripts, or in a dedicated user-defined application created with building blocks from the platform.
2.2. Anatomical models
For simulations involving anatomical models, either the existing virtual population (ViP)  models can be used, or new models can be segmented on the basis of medical image data (e.g. MRI, CT). The ViP models consist of (currently) 10 detailed anatomical human models based on MRI data from healthy volunteers (figure 1). These surface models distinguish between 80 and more than 200 different organs and tissue entities and have, to date, been used in simulations by over 300 groups worldwide. An extensive literature-based tissue parameter database—dielectric properties (dispersion and low frequency), thermal properties, perfusion and density information—is available online and provides recommended parameter values as well as information about the spread of reported measurement values. This database will gradually include more physical and physiological tissue parameters (e.g. mechanical and acoustic properties). An associated online forum allows users to discuss these properties and dynamically extend the database. In addition to these basic human anatomical models, models of pregnant women at different gestational stages (three/seven/nine months) were generated by incorporating available image data from foetuses into a modified version of the female adult model. Various non-human models are available as well (pig, dog, rats, mice).
For segmentation, a dedicated software that allows coupling highly automated segmentation methods (e.g. multi-modal k-means and gamma-method-based clustering, competitive fuzzy connectedness) with highly interactive ones (e.g. intelligent scissors, interactive watershed transformation, brush) was developed. The toolbox also offers filters for noise reduction, speckle removal, label field smoothing, topologically flexible interpolation, skin adding, etc. The resulting label fields can be either transformed into triangulated surfaces (including volume-preserving smoothing, feature-preserving simplification and mesh regularization) or directly used to generate optimized, unstructured tetrahedral meshes with topologically conforming tissue interfaces.
2.4. Posing and morphing
Posing and morphing functionalities are available to parametrize the posture and shape of the anatomical models while preserving realistic internal organ placement and tissue distributions. To pose the models in different postures, two approaches were implemented. The first method makes use of a volume-preserving, skeleton- and influence-region-based approach from the field of computer graphics. This method allows articulating the models and real-time posing. The second approach uses a mechanical simulation (see, §3.4 for more details) with actively prescribed bone movement and passively deformable soft tissue to obtain more realistic joint geometries. A similar approach can be used to further parametrize the models with regard to BMI and thus extend the population coverage. Finite element method (FEM) simulations are performed in which the bones are treated as rigid, the fatty tissues are assigned a growing or shrinking force, and the remaining soft tissues deform passively. In addition to this physics-based morphing approach, an interactive interpolation-based approach was developed, where the user can move control points of various widgets to interactively define a deformation field and change, for example, the shape of internal organs locally. This method can be used to mimic breathing, which could also be simulated using the FEM-based approach, provided the rib cage is suitably articulated.
For simulations involving implanted devices (e.g. pacemakers), either the CAD files can be imported or geometrical modelling can be used together with the anatomical models.
Routines that allow the discretization of the models for various numerical methods were developed. Those routines include homogeneous voxeling, voxeling on arbitrary rectilinear meshes (including conformal subcell generation) with robust intersection testing, and the generation of high-quality, body-fitted, unstructured tetrahedral meshes from CAD or segmentation label fields. For the latter, different methods are available: an advancing front (with optimization), a cut-cell octree (including smoothing) and a Delaunay-based optimization method. Various smoothing and refinement strategies can be used.
A multitude of solvers are available for (i) electromagnetism (harmonic and transient finite-differences time-domain method (FDTD), various quasi-static and low-frequency solvers as well as mode-matching; support for thin sheets, nonlinear materials, circuit co-simulation etc.), (ii) flow (fully implicit FEM and pressure-correction finite volume method; incompressible flow only), (iii) acoustics (full-wave FDTD, nonlinear Westervelt–Lighthill equation), (iv) temperature (FDTD and steady-state; discrete vasculature, thermoregulation, phase changes including latent heat, thin structure models, etc.), (v) CRD (FEM), (vi) solid mechanics (FEM), (vii) pipeline networks (for vasculature mostly), and (viii) neuron dynamics (including interface to NEURON ). All solvers are optimized for large simulations that include complex, inhomogeneous tissue distributions and contain special physical/biological models of particular relevance to simulations in the human body (e.g. multi-scale perfusion and thermoregulation models). Currently, only weak (iterative) coupling of solvers is supported, but a strongly coupled fluid–structure interaction solver is in development. An extended coupling framework with scheduling, data exchange (including interpolation and expression parsing) and load balancing functionality (e.g. on remotely accessible computational resources) is being developed.
2.8. High-performance computing
Most solvers support HPC on various platforms. This includes hardware acceleration (single and multi-GPU) as well as parallelization (message passing interface (MPI), OpenMP). The OpenMP solvers scale well on multi-core, shared memory machines (60–80% parallel efficiency depending on application; tested on up to 12 cores). The MPI parallelization of both the unstructured mesh solvers and the implicit structured mesh solvers is enabled by PETSc . This provides cross-platform support for a large range of hardware architectures, ranging from multi-core desktop PCs to clusters and supercomputers. The flow simulations described below in §3.3 scale-up to several hundred nodes before the MPI version becomes limited by communication. Beside the multiple preconditioners already offered by PETSc, additional preconditioners (e.g. Schur complement, algebraic multi-grid) are available. A generalized framework has been developed that can be used to rapidly implement new HPC-enhanced FEM solvers by simply overloading the coefficient calculation routines. GPU-based hardware acceleration is used for transient simulations on structured, rectilinear meshes. Currently, the use of hybrid MPI–OpenMP approaches is investigated in the context of a general linear solver for indefinite systems.
2.9. Post-processing, optimization and analysis
The simulation platform is complemented by post-processing tools for visualization and analysis, and by an optimization and sensitivity analysis framework.
The main focus areas are design and optimization of medical devices, basic research (investigation of mechanisms/interactions), safety assessment and (patient-specific) treatment planning. The simulation platform has been applied to a large number of relevant problems. Four examples are described below.
3.1. Patient-specific treatment planning and applicator design
To illustrate how patient-specific modelling, treatment planning and applicator design can be achieved with the simulation platform, an example from hyperthermia cancer treatment is presented. Hyperthermia is a cancer treatment modality that aims at selectively heating the tumour region, resulting in preferential killing of tumour cells . For deep-seated tumour regions, RF energy from an antenna array is generally focused into the tumour to heat it. Treatment planning is required to insure sufficient tumour coverage—surviving tumour cells are likely to result in recurrences—and to avoid hot-spots in healthy tissue. Such hot-spots might cause pain or damage in sensitive regions and be treatment limiting. The complex energy deposition (affected by the large dielectric inhomogeneity of the human body and by anatomical structures with dimensions of the order of magnitude of the wavelength) together with the strong influence of local perfusion on the temperature distribution make treatment planning a necessity. In addition to treatment planning, modelling can be used for the development and optimization of novel applicators, for quality assurance and educational purposes, and safety assessment.
A comprehensive hyperthermia treatment planning platform has been implemented using the simulation framework and is currently used in clinics . Patient image data can be segmented by means of the functionality described in §2.3, and a Python-script-based wizard helps to interactively position the patient model in a phased antenna array applicator model. The wizard also takes care of tissue/material properties assignment, discretization, simulation, optimization and report generation. The hardware accelerated (GPU) EM FDTD solver is used to perform EM simulations for each antenna. Various general (e.g. genetic algorithms) as well as dedicated (e.g. generalized eigenvector method) optimizers are then used to find the optimal steering parameters (figure 2). The transient heating that results from the deposited energy is calculated, and thermoregulation, local perfusion effects and the impact of major vessels can be taken into consideration. To assess the expected treatment outcome, the transient tissue temperature at each location can be translated into tissue damage based on an Arrhenius model. Alternatively it can be used to predict the administered thermal dose according to the cumulative equivalent minutes at 43°C (CEM43, ) concept. Real-time re-optimization based on patient feedback has been achieved with a Pareto multi-objective optimization approach. When the patient complains about hot-spots or when measurements indicate insufficient heating of some parts of the tumour, the different objectives (heating of various tumour regions, avoiding exposure of sensitive healthy tissue) can be reweighted to modify treatment parameters in real-time therby bringing treatment planning into the treatment room.
The platform has also been used to develop a novel head and neck applicator [10,11] and to assess the safety of hyperthermia in the presence of metallic implants. It was validated on the basis of dosimetric measurements and clinical information. Further experimental validation, especially of the perfusion model, is planned, and the segmentation process must be accelerated to facilitate routine clinical use of hyperthermia treatment planning. The large number of tissues that need to be distinguished in order to handle the widely varying dielectric and perfusion parameters render the latter complicated. Using image-based tissue parameter estimation supported by the simulation platform could reduce the segmentation requirements.
Similar methodology has been used to assess the safety of MR exposure in the presence and absence of implants to provide relevant information for revising MR safety standards, for implant safety certification, and for the investigation of novel parallel-transmit coil technology.
3.2. Tissue modelling and basic research
To illustrate the use of the simulation platform for basic research, mechanism investigation and tissue modelling, an example from the field of tumour growth modelling is presented. The interplay of factors at the tissue, cellular and subcellular levels has been considered in simulations of solid tumour growth . The model takes into account neoplastic tissue growth, blood and oxygen transport (by means of a hybrid approach that combines continuum scale reaction–diffusion equations and microscale transport in a developing pipeline vessel network) and angiogenic sprouting. Angiogenesis is treated by explicitly modelling sprouting vessel tips guided by chemotaxis along angiogenic growth factor (AGF) gradients. AGF is released in hypoxic regions of the tumour. Tumour and endothelial cell proliferation, migration, quiescence and apoptosis depend on the presence of nutrients as well as on local mechanical stress. As a result of tissue growth and shrinkage, mechanical forces are determined that cause geometry deformation.
The numerical implementation involves a coupled system of (i) CRD equations and (ii) mechanical simulations performed with the FEM solvers as well as (iii) a pipeline solver (for momentum and mass transport as well as oxygen delivery through the vessel wall). The coupling has to consider that the various processes occur on different temporal scales ranging from nanoseconds (molecular diffusion) to days (cell proliferation). Therefore, some processes are solved in a quasi-static manner whereas others are treated as transient.
The model is capable of correctly reproducing some major solid tumour features such as the necrotic core, the brush-shaped vasculature capsule (figure 3), or the incident angle of exophytic tumours (provided suitable nonlinear tissue elasticity models are used). It has been used to study the impact of various treatment approaches such as vessel embolization, the administration of angiogenesis inhibitors such as angiostatin, or hyperthermia (as a single treatment modality or in combination with chemotherapy). When modelling hyperthermia treatment of cancer, the tumour model was extended by considering the temperature dependence of the perfusion (and related nutrient delivery), the direct temperature-induced apoptosis as well as indirect cell killing owing to the enhanced activity of a chemotherapeutic agent such as doxorubicin . The model was also used to assess whether the local administration of angiogenesis inhibitors can lead to effective tumour treatment or only causes the vasculature to circumvent the treated area .
In a next step, it will be important to validate the model beyond the phenomenological, macroscopic level and to continue including other relevant aspects such as interstitial fluid pressure to increase the realism of the model.
3.3. Image-based modelling
As an example of image-based modelling, an investigation of the magneto-haemodynamic (MHD) effect as a biomarker is presented. Here image data were not only used to provide geometry, but as well boundary condition information. Blood flow inside a strong magnetic field (e.g. an MRI magnet) leads to surface potential changes that are measurable as changes in the electrocardiogram (ECG). These MHD effect-related changes can be problematic in the context of ECG-based triggering. To investigate whether such ECG changes can be used as biomarkers for blood flow features (e.g. to measure cardiac output or recirculation), simulations based on medical image data from a volunteer have been performed and compared with experimental measurements .
Both the geometry and the transient inflow boundary conditions have been extracted from MRI measurements and used to perform flow simulations of the aorta and vena cava. The simulations used the FEM solver with an approximate Schur-complement preconditioner for incompressible flow (figure 4). For the outflow boundary conditions, developed flow and image-based volume rates have been imposed. The geometry was treated as rigid, neglecting wall pulsation. The obtained transient flow results were validated by comparisons with both MR flow measurements in selected cross-sections and four-dimensional measurement data of the aortic arc region. The blood flow-induced surface potential change could then be determined by means of the low-frequency EM solver (figure 4) with a special velocity-related source term.
ECG measurements have been performed both inside and outside the MRI magnet at 11 locations on the frontal torso of the volunteer. By averaging the ECG signal over approximately 100 heart beats (after aligning the R peak of the QRS complex) and subtracting the measurements taken outside of the magnet from those taken inside, the MHD effect-related signal could be extracted. This signal was found to be strongly position-dependent with highly variable, but reproducibly measurable characteristics. This indicates the possibility of deducing location specific information on blood flow.
The simulated MHD signal agrees with major characteristics of the measured signal (temporal location of main peak, magnitude, variation across chest and along torso ) except in the vicinity of the heart, where the impact of the non-simulated blood flow in the heart is important. The MHD signal is primarily influenced by the aorta. However, more vessels and better boundary conditions are needed to analyse the finer details of the response. The simulations allow the study of how flow features correlate with the MHD signal and which aspects of the MHD signal are caused by a specific subvolume of the flow (e.g. the flow in the aortic arc, at a bifurcation or in a vortex zone).
3.4. Anatomical models and morphing
To illustrate the anatomical models and their parametrization, physics-based morphing of the anatomical models is described. Mechanical simulations have been performed to parameterize the ViP anatomical models with regard to fat content, to increase the population coverage: the bones are treated as rigid, while the fatty tissue is assigned a growing or shrinking force. The remaining soft tissues are deformed passively (figure 1) on the basis of realistic mechanical properties. This approach is motivated by the need for physiological fidelity and requires mechanical simulations involving over 100 000 000 elements that can be efficiently executed on a desktop PC with 144 Gb of RAM by the HPC-enabled FEM solver. During the mechanical simulations, the appearance of self-intersections is carefully monitored and remeshing is applied if necessary. An efficient linear-elastic, small-strain-approximation mechanical solver is used for small deformations, and a Neohookean tissue model is typically used for larger ones. More complex nonlinear material models based on generic energy functions are currently being integrated.
These parametrized versions extend the range of the population covered by the ViP models and are valuable, for example, when assessing exposure safety of the general public. The same technique is also used to realistically change the posture of the models by repositioning their bones via pivots in the joints and passively deforming the soft tissue by means of a mechanical simulation. Compared with the aforementioned influence-region-based approach (§3.4) this offers improved realism especially in articulations and when mimicking large posture changes. A similar approach to the fat morphing could be used to parametrize changes in muscle mass, possibly using anisotropic growth forces based on fibre orientation. The morphing approach has been used in the context of MRI RF-safety assessment to generate a multitude of typical scanning postures from image data taken from a volunteer in supine position. When generating a prone position, it can be necessary to include gravitational forces as well as contact (MRI table) in the simulation.
A multi-physics simulation platform optimized for CLS has been developed, with device development and optimization, safety assessment, basic research, and treatment planning as primary application areas. It strongly supports medical image data-based modelling and simulations involving complex anatomical models. HPC-enabled solvers enable simulations of large, realistic problems on a wide range of computational architectures. Various solvers (e.g. neuronal, thermal, ultrasound) have been specially developed for simulations in living tissue. The simulation platform integrates VTK visualization and Python scripting and allows users to create their own applications and plug-ins. It has been used to tackle a large range of relevant problems such as those presented in this paper, i.e. hyperthermia cancer treatment planning, tumour growth modelling, studying the use of the MHD effect as a biomarker, and physics-based morphing of anatomical models.
It is hoped that the availability of such simulation platforms will help modelling become a more widely used and beneficial tool in the medical technology and life sciences fields, thereby increasing the understanding, accelerating the development of treatments and devices, and enhancing safety for the benefit of the individual patients.
Various components of this work have been supported by different organizations: CTI, SNSF CO-ME, SciEx, ISJRP. Various partners have contributed to the development: ETHZ, Erasmus MC, FDA, MIT, EPFL, UniBas, University Hospital Zurich, CABMM, CSCS, UZH, SPEAG, ZMT.
One contribution of 25 to a Theme Issue ‘The virtual physiological human: integrative approaches to computational biomedicine’.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.