Cardiac digital twins are computational tools capturing key functional and anatomical characteristics of patient hearts for investigating disease phenotypes and predicting responses to therapy. When paired with large-scale computational resources and large clinical datasets, digital twin technology can enable virtual clinical trials on virtual cohorts to fast-track therapy development. Here, we present an open-source automated pipeline for personalising ventricular electrophysiological function based on routinely acquired magnetic resonance imaging (MRI) data and the standard 12-lead electrocardiogram (ECG). Using MRI-based anatomical models, a sequential Monte-Carlo approximate Bayesian computational inference method is extended to infer electrical activation and repolarisation characteristics from the ECG. Fast simulations are conducted with a reaction-Eikonal model, including the Purkinje network and biophysically-detailed subcellular ionic current dynamics for repolarisation. For each patient, parameter uncertainty is represented by inferring an envelope of plausible ventricular models rather than a single one, which means that parameter uncertainty can be propagated to therapy evaluation. Furthermore, we have developed techniques for translating from reaction-Eikonal to monodomain simulations, which allows more realistic simulations of cardiac electrophysiology. The pipeline is demonstrated in three healthy subjects, where our inferred pseudo-diffusion reaction-Eikonal models reproduced the patient's ECG with a median Pearson's correlation coefficient of 0.9, and then translated to monodomain simulations with a median correlation coefficient of 0.84 across all subjects. We then demonstrate our digital twins for virtual evaluation of Dofetilide with uncertainty quantification. These evaluations using our cardiac digital twins reproduced dose-dependent QTc and T peak to T end prolongations that are in keeping with large population drug response data. The methodologies for cardiac digital twinning presented here are a step towards personalised virtual therapy testing and can be scaled to generate virtual populations for clinical trials to fast-track therapy evaluation. The tools developed for this paper are open-source, documented, and made publicly available.

Multiple myeloma (MM) is a plasma cell cancer that occurs in the bone marrow. A leading treatment for MM is the monoclonal antibody Daratumumab, targeting the CD38 receptor, which is highly overexpressed in myeloma cells. In this work we model drug resistance via loss of CD38 expression, which is a proposed mechanism of resistance to Daratumumab treatment. We develop an ODE model that includes drug resistance via two mechanisms: a direct effect in which CD38 expression is lost without cell death in response to Daratumumab, and an indirect effect in which CD38 expression switches on and off in the cancer cells; myeloma cells that do not express CD38 have lower fitness but are shielded from the drug action. The model also incorporates competition with healthy cells, death of healthy cells due to off-target drug effects, and a Michaelis-Menten type immune response. Using optimal control theory, we study the effect of the drug resistance mechanisms and the off-target drug effect on the optimal treatment regime. We identify a general increase in the duration and costs of optimal treatment, as a result of these added mechanisms. Several distinct optimal treatment regimes are identified within the parameter space.

Cardiac in silico clinical trials can virtually assess the safety and efficacy of therapies using human-based modelling and simulation. These technologies can provide mechanistic explanations for clinically observed pathological behaviour. Designing virtual cohorts for in silico trials requires exploiting clinical data to capture the physiological variability in the human population. The clinical characterisation of ventricular activation and the Purkinje network is challenging, especially non-invasively. Our study aims to present a novel digital twinning pipeline that can efficiently generate and integrate Purkinje networks into human multiscale biventricular models based on subject-specific clinical 12-lead electrocardiogram and magnetic resonance recordings. Essential novel features of the pipeline are the human-based Purkinje network generation method, personalisation considering ECG R wave progression as well as QRS morphology, and translation from reduced-order Eikonal models to equivalent biophysically-detailed monodomain ones. We demonstrate ECG simulations in line with clinical data with clinical image-based multiscale models with Purkinje in four control subjects and two hypertrophic cardiomyopathy patients (simulated and clinical QRS complexes with Pearson's correlation coefficients > 0.7). Our methods also considered possible differences in the density of Purkinje myocardial junctions in the Eikonal-based inference as regional conduction velocities. These differences translated into regional coupling effects between Purkinje and myocardial models in the monodomain formulation. In summary, we demonstrate a digital twin pipeline enabling simulations yielding clinically consistent ECGs with clinical CMR image-based biventricular multiscale models, including personalised Purkinje in healthy and cardiac disease conditions.

Fibrosis, a pathological increase in extracellular matrix proteins, is a significant health issue that hinders the function of many organs in the body, in some cases fatally. In the heart, fibrosis impacts on electrical propagation in a complex and poorly predictable fashion, potentially serving as a substrate for dangerous arrhythmias. Individual risk depends on the spatial manifestation of fibrotic tissue, and learning the spatial arrangement on the fine scale in order to predict these impacts still relies upon invasive ex vivo procedures. As a result, the effects of spatial variability on the symptomatic impact of cardiac fibrosis remain poorly understood. In this work, we address the issue of availability of such imaging data via a computational methodology for generating new realisations of cardiac fibrosis microstructure. Using the Perlin noise technique from computer graphics, together with an automated calibration process that requires only a single training image, we demonstrate successful capture of collagen texturing in four types of fibrosis microstructure observed in histological sections. We then use this generator to quantitatively analyse the conductive properties of these different types of cardiac fibrosis, as well as produce three-dimensional realisations of histologically-observed patterning. Owing to the generator's flexibility and automated calibration process, we also anticipate that it might be useful in producing additional realisations of other physiological structures.

The realisation of precision cardiology requires novel techniques for the non-invasive characterisation of individual patients’ cardiac function to inform therapeutic and diagnostic decision-making. Both electrocardiography and imaging are used for the clinical diagnosis of cardiac disease. The integration of multi-modal datasets through advanced computational methods could enable the development of the cardiac ‘digital twin’, a comprehensive virtual tool that mechanistically reveals a patient's heart condition from clinical data and simulates treatment outcomes. The adoption of cardiac digital twins requires the non-invasive efficient personalisation of the electrophysiological properties in cardiac models. This study develops new computational techniques to estimate key ventricular activation properties for individual subjects by exploiting the synergy between non-invasive electrocardiography, cardiac magnetic resonance (CMR) imaging and modelling and simulation. More precisely, we present an efficient sequential Monte Carlo approximate Bayesian computation-based inference method, integrated with Eikonal simulations and torso-biventricular models constructed based on clinical CMR imaging. The method also includes a novel strategy to treat combined continuous (conduction speeds) and discrete (earliest activation sites) parameter spaces and an efficient dynamic time warping-based ECG comparison algorithm. We demonstrate results from our inference method on a cohort of twenty virtual subjects with cardiac ventricular myocardial-mass volumes ranging from 74 cm3 to 171 cm3 and considering low versus high resolution for the endocardial discretisation (which determines possible locations of the earliest activation sites). Results show that our method can successfully infer the ventricular activation properties in sinus rhythm from non-invasive epicardial activation time maps and ECG recordings, achieving higher accuracy for the endocardial speed and sheet (transmural) speed than for the fibre or sheet-normal directed speeds.

The titre of virus in a dengue patient and the duration of this viraemia has a profound effect on whether or not a mosquito will become infected when it feeds on the patient and this, in turn, is a key driver of the magnitude of a dengue outbreak. The assessment of the heterogeneity of viral dynamics in dengue-infected patients and its precise treatment are still uncertain. Infection onset, patient physiology and immune response are thought to play major roles in the development of the viral load. Research has explored the interference and spontaneous generation of defective virus particles, but have not examined both the antibody and defective particles during natural infection. We explore the intrinsic variability in the within-host dynamics of viraemias for a population of patients using the method of population of models (POMs). A dataset from 208 patients is used to initially calibrate 20,000 models for the infection kinetics for each of the four dengue virus serotypes. The calibrated POMs suggests that naturally generated defective particles may interfere with the viraemia, but the generated defective virus particles are not adequate to reduce high fever and viraemia duration. The effect of adding excess defective dengue virus interfering particles to patients as a therapeutic is evaluated using the calibrated POMs in a bang-bang (on-off or two-step) optimal control setting. Bang-bang control is a class of binary feedback control that turns either 'ON' or 'OFF' at different time points, determined by the system feedback. Here, the bang-bang control estimates the mathematically optimal dose and duration of the intervention for each model in the POM set.

Welcome to the 2016 IEEE International Conference on Bioinformatics and Biomedicine (IEEE BIBM 2016) being held in the Shenzhen, China from December 15–18, 2016. On behalf of the IEEE BIBM 2016 Organizing Team, we would like to thank you for your participation and hope you enjoy the conference.

Diffusion in biological tissues is known to be hindered by the structural complexity of the underlying medium. In the heart, improved characterisation on how this complexity influences acquired diffusion weighted signals is key to advancing our interpretation of diffusion magnetic resonance imaging, as well as to propose novel biomarkers to further characterise myocardial microstructure. In this work, we propose stretched Mittag-Leffler signal decay models for the quantification of the anomalous decay observed in acquired diffusion weighted signals. Our results, analysed in ex vivo healthy, fixed rat ventricles, indicate that such a representation suffices to capture the anomalous signal decay observed in the myocardial syncytium. The subdiffusive order of signal decay is shown to encode independent information to that encapsulated by standard diffusion tensor metrics, and thus may provide additional information on tissue microstructure. Moreover, subdiffusion gradients are shown to be indicative of the total structural heterogeneity spanning the left ventricular wall, which includes progressive myolaminae branching and spatially varying densities of perimysial collagen, microvasculature and adipose tissue. The proposed approach may therefore have important implications for the characterisation of tissue microstructure, both in cardiac and other tissue types.

Since we still know very little about stem cells in their natural environment, it is useful to explore their dynamics through modelling and simulation, as well as experimentally. Most models of stem cell systems are based on deterministic differential equations that ignore the natural heterogeneity of stem cell populations. This is not appropriate at the level of individual cells and niches, when randomness is more likely to affect dynamics. In this paper, we introduce a fast stochastic method for simulating a metapopulation of stem cell niche lineages, that is, many sub-populations that together form a heterogeneous metapopulation, over time. By selecting the common limiting timestep, our method ensures that the entire metapopulation is simulated synchronously. This is important, as it allows us to introduce interactions between separate niche lineages, which would otherwise be impossible. We expand our method to enable the coupling of many lineages into niche groups, where differentiated cells are pooled within each niche group. Using this method, we explore the dynamics of the haematopoietic system from a demand control system perspective. We find that coupling together niche lineages allows the organism to regulate blood cell numbers as closely as possible to the homeostatic optimum. Furthermore, coupled lineages respond better than uncoupled ones to random perturbations, here the loss of some myeloid cells. This could imply that it is advantageous for an organism to connect together its niche lineages into groups. Our results suggest that a potential fruitful empirical direction will be to understand how stem cell descendants communicate with the niche and how cancer may arise as a result of a failure of such communication.

It is widely accepted that lipid rafts promote receptor clustering and thereby facilitate signaling transduction. The role of lipid rafts in inducing and promoting receptor accumulation within the cell membrane has been explored by several computational and experimental studies. However, it remains unclear whether lipid rafts influence the recruitment and binding of proteins from the cytosol as well. To provide an answer to this question a spatial membrane model has been developed based on cellular automata. Our results indicate that lipid rafts indeed influence protein receptor bindings. In particular processes with slow dissociation and binding kinetics are promoted by lipid rafts, whereas fast binding processes are slightly hampered. However, the impact depends on a variety of parameters, such as the size and mobility of the lipid rafts, the induced slow down of receptors within rafts, and also the dissociation and binding kinetics of the cytosolic proteins. Thus, for any individual signaling pathway the influence of lipid rafts on protein binding might be different. To facilitate analyzing this influence given a specific pathway, our approach has been generalized into LiRaM, a modeling and simulation tool for lipid rafts models.

Markov clustering (MCL) is becoming a key algorithm within bioinformatics for determining clusters in networks. However,with increasing vast amount of data on biological networks, performance and scalability issues are becoming a critical limiting factor in applications. Meanwhile, GPU computing, which uses CUDA tool for implementing a massively parallel computing environment in the GPU card, is becoming a very powerful, efficient, and low-cost option to achieve substantial performance gains over CPU approaches. The use of on-chip memory on the GPU is efficiently lowering the latency time, thus, circumventing a major issue in other parallel computing environments, such as MPI. We introduce a very fast Markov clustering algorithm using CUDA (CUDA-MCL) to perform parallel sparse matrix-matrix computations and parallel sparse Markov matrix normalizations, which are at the heart of MCL. We utilized ELLPACK-R sparse format to allow the effective and fine-grain massively parallel processing to cope with the sparse nature of interaction networks data sets in bioinformatics applications. As the results show, CUDA-MCL is significantly faster than the original MCL running on CPU. Thus, large-scale parallel computation on off-the-shelf desktop-machines, that were previously only possible on supercomputing architectures, can significantly change the way bioinformaticians and biologists deal with their data.

Chemotaxis plays a crucial role in many biological processes, including nervous system development. However, fundamental physical constraints limit the ability of a small sensing device such as a cell or growth cone to detect an external chemical gradient. One of these is the stochastic nature of receptor binding, leading to a constantly fluctuating binding pattern across the cell's array of receptors. This is analogous to the uncertainty in sensory information often encountered by the brain at the systems level. Here we derive analytically the Bayes-optimal strategy for combining information from a spatial array of receptors in both one and two dimensions to determine gradient direction. We also show how information from more than one receptor species can be optimally integrated, derive the gradient shapes that are optimal for guiding cells or growth cones over the longest possible distances, and illustrate that polarized cell behavior might arise as an adaptation to slowly varying environments. Together our results provide closed-form predictions for variations in chemotactic performance over a wide range of gradient conditions.

One of the most important aspects of Computational Cell Biology is the understanding of the complicated dynamical processes that take place on plasma membranes. These processes are often so complicated that purely temporal models cannot always adequately capture the dynamics. On the other hand, spatial models can have large computational overheads. In this article, we review some of these issues with respect to chemistry, membrane microdomains and anomalous diffusion and discuss how to select appropriate modelling and simulation paradigms based on some or all the following aspects: discrete, continuous, stochastic, delayed and complex spatial processes.

MOTIVATION: Disulfide bonds are primary covalent crosslinks between two cysteine residues in proteins that play critical roles in stabilizing the protein structures and are commonly found in extracy-toplasmatic or secreted proteins. In protein folding prediction, the localization of disulfide bonds can greatly reduce the search in conformational space. Therefore, there is a great need to develop computational methods capable of accurately predicting disulfide connectivity patterns in proteins that could have potentially important applications. RESULTS: We have developed a novel method to predict disulfide connectivity patterns from protein primary sequence, using a support vector regression (SVR) approach based on multiple sequence feature vectors and predicted secondary structure by the PSIPRED program. The results indicate that our method could achieve a prediction accuracy of 74.4% and 77.9%, respectively, when averaged on proteins with two to five disulfide bridges using 4-fold cross-validation, measured on the protein and cysteine pair on a well-defined non-homologous dataset. We assessed the effects of different sequence encoding schemes on the prediction performance of disulfide connectivity. It has been shown that the sequence encoding scheme based on multiple sequence feature vectors coupled with predicted secondary structure can significantly improve the prediction accuracy, thus enabling our method to outperform most of other currently available predictors. Our work provides a complementary approach to the current algorithms that should be useful in computationally assigning disulfide connectivity patterns and helps in the annotation of protein sequences generated by large-scale whole-genome projects. AVAILABILITY: The prediction web server and Supplementary Material are accessible at http://foo.maths.uq.edu.au/~huber/disulfide

MOTIVATION: Kinetic rate in gene expression is a key measurement of the stability of gene products and gives important information for the reconstruction of genetic regulatory networks. Recent developments in experimental technologies have made it possible to measure the numbers of transcripts and protein molecules in single cells. Although estimation methods based on deterministic models have been proposed aimed at evaluating kinetic rates from experimental observations, these methods cannot tackle noise in gene expression that may arise from discrete processes of gene expression, small numbers of mRNA transcript, fluctuations in the activity of transcriptional factors and variability in the experimental environment. RESULTS: In this paper, we develop effective methods for estimating kinetic rates in genetic regulatory networks. The simulated maximum likelihood method is used to evaluate parameters in stochastic models described by either stochastic differential equations or discrete biochemical reactions. Different types of non-parametric density functions are used to measure the transitional probability of experimental observations. For stochastic models described by biochemical reactions, we propose to use the simulated frequency distribution to evaluate the transitional density based on the discrete nature of stochastic simulations. The genetic optimization algorithm is used as an efficient tool to search for optimal reaction rates. Numerical results indicate that the proposed methods can give robust estimations of kinetic rates with good accuracy.

Systems biology is based on computational modelling and simulation of large networks of interacting components. Models may be intended to capture processes, mechanisms, components and interactions at different levels of fidelity. Input data are often large and geographically disperse, and may require the computation to be moved to the data, not vice versa. In addition, complex system-level problems require collaboration across institutions and disciplines. Grid computing can offer robust, scaleable solutions for distributed data, compute and expertise. We illustrate some of the range of computational and data requirements in systems biology with three case studies: one requiring large computation but small data (orthologue mapping in comparative genomics), a second involving complex terabyte data (the Visible Cell project) and a third that is both computationally and data-intensive (simulations at multiple temporal and spatial scales). Authentication, authorisation and audit systems are currently not well scalable and may present bottlenecks for distributed collaboration particularly where outcomes may be commercialised. Challenges remain in providing lightweight standards to facilitate the penetration of robust, scalable grid-type computing into diverse user communities to meet the evolving demands of systems biology.

BACKGROUND: The residue-wise contact order (RWCO) describes the sequence separations between the residues of interest and its contacting residues in a protein sequence. It is a new kind of one-dimensional protein structure that represents the extent of long-range contacts and is considered as a generalization of contact order. Together with secondary structure, accessible surface area, the B factor, and contact number, RWCO provides comprehensive and indispensable important information to reconstructing the protein three-dimensional structure from a set of one-dimensional structural properties. Accurately predicting RWCO values could have many important applications in protein three-dimensional structure prediction and protein folding rate prediction, and give deep insights into protein sequence-structure relationships. RESULTS: We developed a novel approach to predict residue-wise contact order values in proteins based on support vector regression (SVR), starting from primary amino acid sequences. We explored seven different sequence encoding schemes to examine their effects on the prediction performance, including local sequence in the form of PSI-BLAST profiles, local sequence plus amino acid composition, local sequence plus molecular weight, local sequence plus secondary structure predicted by PSIPRED, local sequence plus molecular weight and amino acid composition, local sequence plus molecular weight and predicted secondary structure, and local sequence plus molecular weight, amino acid composition and predicted secondary structure. When using local sequences with multiple sequence alignments in the form of PSI-BLAST profiles, we could predict the RWCO distribution with a Pearson correlation coefficient (CC) between the predicted and observed RWCO values of 0.55, and root mean square error (RMSE) of 0.82, based on a well-defined dataset with 680 protein sequences. Moreover, by incorporating global features such as molecular weight and amino acid composition we could further improve the prediction performance with the CC to 0.57 and an RMSE of 0.79. In addition, combining the predicted secondary structure by PSIPRED was found to significantly improve the prediction performance and could yield the best prediction accuracy with a CC of 0.60 and RMSE of 0.78, which provided at least comparable performance compared with the other existing methods. CONCLUSION: The SVR method shows a prediction performance competitive with or at least comparable to the previously developed linear regression-based methods for predicting RWCO values. In contrast to support vector classification (SVC), SVR is very good at estimating the raw value profiles of the samples. The successful application of the SVR approach in this study reinforces the fact that support vector regression is a powerful tool in extracting the protein sequence-structure relationship and in estimating the protein structural profiles from amino acid sequences.

BACKGROUND: The majority of peptide bonds in proteins are found to occur in the trans conformation. However, for proline residues, a considerable fraction of Prolyl peptide bonds adopt the cis form. Proline cis/trans isomerization is known to play a critical role in protein folding, splicing, cell signaling and transmembrane active transport. Accurate prediction of proline cis/trans isomerization in proteins would have many important applications towards the understanding of protein structure and function. RESULTS: In this paper, we propose a new approach to predict the proline cis/trans isomerization in proteins using support vector machine (SVM). The preliminary results indicated that using Radial Basis Function (RBF) kernels could lead to better prediction performance than that of polynomial and linear kernel functions. We used single sequence information of different local window sizes, amino acid compositions of different local sequences, multiple sequence alignment obtained from PSI-BLAST and the secondary structure information predicted by PSIPRED. We explored these different sequence encoding schemes in order to investigate their effects on the prediction performance. The training and testing of this approach was performed on a newly enlarged dataset of 2424 non-homologous proteins determined by X-Ray diffraction method using 5-fold cross-validation. Selecting the window size 11 provided the best performance for determining the proline cis/trans isomerization based on the single amino acid sequence. It was found that using multiple sequence alignments in the form of PSI-BLAST profiles could significantly improve the prediction performance, the prediction accuracy increased from 62.8% with single sequence to 69.8% and Matthews Correlation Coefficient (MCC) improved from 0.26 with single local sequence to 0.40. Furthermore, if coupled with the predicted secondary structure information by PSIPRED, our method yielded a prediction accuracy of 71.5% and MCC of 0.43, 9% and 0.17 higher than the accuracy achieved based on the singe sequence information, respectively. CONCLUSION: A new method has been developed to predict the proline cis/trans isomerization in proteins based on support vector machine, which used the single amino acid sequence with different local window sizes, the amino acid compositions of local sequence flanking centered proline residues, the position-specific scoring matrices (PSSMs) extracted by PSI-BLAST and the predicted secondary structures generated by PSIPRED. The successful application of SVM approach in this study reinforced that SVM is a powerful tool in predicting proline cis/trans isomerization in proteins and biological sequence analysis.

Discrete stochastic simulations are a powerful tool for understanding the dynamics of chemical kinetics when there are small-to-moderate numbers of certain molecular species. In this paper we introduce delays into the stochastic simulation algorithm, thus mimicking delays associated with transcription and translation. We then show that this process may well explain more faithfully than continuous deterministic models the observed sustained oscillations in expression levels of hes1 mRNA and Hes1 protein.

Recent advances in gene-expression profiling technologies provide large amounts of gene expression data. This raises the possibility for a functional understanding of genome dynamics by means of mathematical modelling. As gene expression involves intrinsic noise, stochastic models are essential for better descriptions of gene regulatory networks. However, stochastic modelling for large scale gene expression data sets is still in the very early developmental stage. In this paper we present some stochastic models by introducing stochastic processes into neural network models that can describe intermediate regulation for large scale gene networks. Poisson random variables are used to represent chance events in the processes of synthesis and degradation. For expression data with normalized concentrations, exponential or normal random variables are used to realize fluctuations. Using a network with three genes, we show how to use stochastic simulations for studying robustness and stability properties of gene expression patterns under the influence of noise, and how to use stochastic models to predict statistical distributions of expression levels in population of cells. The discussion suggest that stochastic neural network models can give better description of gene regulatory networks and provide criteria for measuring the reasonableness o mathematical models.

The task of fitting smoothing spline surfaces to meteorological data such as temperature or rainfall observations is computationally intensive. The Generalised Cross Validation (GCV) smoothing algorithm is O(n3) computationally, and memory requirements are 0(n2). Fitting a spline to a moderately sized data set of, for example. 1080 observations and calculating an output surface grid of dimension 220 × 220 involves approximately 5 billion floating point operations, and takes approximately 19 minutes of execution time on a Sun SPARC2 workstation. Since fitting a surface to data collected from the whole of Australia could conceivably involve data sets with approximately 10000 points, and because it is desirable to be able to fit surfaces of at least 1000 data points in 1 to 5 seconds for use in interactive visualisations, it is crucial to be able to take advantage of supercomputing resources. This paper describes the adaptation of the surface fitting program to different supercomputing platforms, and the results achieved.