On the surprising effectiveness of a simple matrix exponential derivative approximation, with application to global SARS-CoV-2.

The continuous-time Markov chain (CTMC) is the mathematical workhorse of evolutionary biology. Learning CTMC model parameters using modern, gradient-based methods requires the derivative of the matrix exponential evaluated at the CTMC's infinitesimal generator (rate) matrix. Motivated by the derivative's extreme computational complexity as a function of state space cardinality, recent work demonstrates the surprising effectiveness of a naive, first-order approximation for a host of problems in computational biology. In response to this empirical success, we obtain rigorous deterministic and probabilistic bounds for the error accrued by the naive approximation and establish a "blessing of dimensionality" result that is universal for a large class of rate matrices with random entries. Finally, we apply the first-order approximation within surrogate-trajectory Hamiltonian Monte Carlo for the analysis of the early spread of Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) across 44 geographic regions that comprise a state space of unprecedented dimensionality for unstructured (flexible) CTMC models within evolutionary biology.

A uniformisation-driven algorithm for inference-related estimation of a phase-type ageing model.

We develop an efficient algorithm to compute the likelihood of the phase-type ageing model. The proposed algorithm uses the uniformisation method to stabilise the numerical calculation. It also utilises a vectorised formula to only calculate the necessary elements of the probability distribution. Our algorithm, with an error's upper bound, could be adjusted easily to tackle the likelihood calculation of the Coxian models. Furthermore, we compare the speed and the accuracy of the proposed algorithm with those of the traditional method using the matrix exponential. Our algorithm is faster and more accurate than the traditional method in calculating the likelihood. Based on our experiments, we recommend using 20 sets of randomly-generated initial values for the optimisation to get a reliable estimate for which the evaluated likelihood is close to the maximum likelihood.

Direct statistical inference for finite Markov jump processes via the matrix exponential.

Given noisy, partial observations of a time-homogeneous, finite-statespace Markov chain, conceptually simple, direct statistical inference is available, in theory, via its rate matrix, or infinitesimal generator, Q , since exp ( Q t ) is the transition matrix over time t. However, perhaps because of inadequate tools for matrix exponentiation in programming languages commonly used amongst statisticians or a belief that the necessary calculations are prohibitively expensive, statistical inference for continuous-time Markov chains with a large but finite state space is typically conducted via particle MCMC or other relatively complex inference schemes. When, as in many applications Q arises from a reaction network, it is usually sparse. We describe variations on known algorithms which allow fast, robust and accurate evaluation of the product of a non-negative vector with the exponential of a large, sparse rate matrix. Our implementation uses relatively recently developed, efficient, linear algebra tools that take advantage of such sparsity. We demonstrate the straightforward statistical application of the key algorithm on a model for the mixing of two alleles in a population and on the Susceptible-Infectious-Removed epidemic model.

Computable upper error bounds for Krylov approximations to matrix exponentials and associated φ -functions.

An a posteriori estimate for the error of a standard Krylov approximation to the matrix exponential is derived. The estimate is based on the defect (residual) of the Krylov approximation and is proven to constitute a rigorous upper bound on the error, in contrast to existing asymptotical approximations. It can be computed economically in the underlying Krylov space. In view of time-stepping applications, assuming that the given matrix is scaled by a time step, it is shown that the bound is asymptotically correct (with an order related to the dimension of the Krylov space) for the time step tending to zero. This means that the deviation of the error estimate from the true error tends to zero faster than the error itself. Furthermore, this result is extended to Krylov approximations of φ -functions and to improved versions of such approximations. The accuracy of the derived bounds is demonstrated by examples and compared with different variants known from the literature, which are also investigated more closely. Alternative error bounds are tested on examples, in particular a version based on the concept of effective order. For the case where the matrix exponential is used in time integration algorithms, a step size selection strategy is proposed and illustrated by experiments.

Exploring inductive linearization for pharmacokinetic-pharmacodynamic systems of nonlinear ordinary differential equations.

Pharmacokinetic-pharmacodynamic systems are often expressed with nonlinear ordinary differential equations (ODEs). While there are numerous methods to solve such ODEs these methods generally rely on time-stepping solutions (e.g. Runge-Kutta) which need to be matched to the characteristics of the problem at hand. The primary aim of this study was to explore the performance of an inductive approximation which iteratively converts nonlinear ODEs to linear time-varying systems which can then be solved algebraically or numerically. The inductive approximation is applied to three examples, a simple nonlinear pharmacokinetic model with Michaelis-Menten elimination (E1), an integrated glucose-insulin model and an HIV viral load model with recursive feedback systems (E2 and E3, respectively). The secondary aim of this study was to explore the potential advantages of analytically solving linearized ODEs with two examples, again E3 with stiff differential equations and a turnover model of luteinizing hormone with a surge function (E4). The inductive linearization coupled with a matrix exponential solution provided accurate predictions for all examples with comparable solution time to the matched time-stepping solutions for nonlinear ODEs. The time-stepping solutions however did not perform well for E4, particularly when the surge was approximated by a square wave. In circumstances when either a linear ODE is particularly desirable or the uncertainty in matching the integrator to the ODE system is of potential risk, then the inductive approximation method coupled with an analytical integration method would be an appropriate alternative.

Dual core processors: Coupled queues: Transient performance evaluation.

High-Performance Computing (HiPC) systems routinely employ multi/many - core processors. Specifically, dual - core processors find many applications in pervasive computing devices. Dual-core processors employ buffers for queueing the incoming jobs. Traditionally, the queues at the processors are assumed to be independent and the queueing system is analyzed in equilibrium for tractability purposes. Queues are modeled using Continuous Time Markov Chains (CTMC's) and the equilibrium performance measures are determined to analyze as well as design the queueing systems. In most interesting cases, the incoming jobs are routed to the queues using the Join the Shortest Queue (JSQ) policy. Thus, with such an adaptive routing algorithm, the two queues are evidently coupled and are not statistically independent. Hence traditional equilibrium performance evaluation doesn't provide realistic performance measures. In this research paper, the two queues associated with buffers in dual-core processors are considered to be coupled. The Coupled Queues are modeled using a Quasi - Birth - and - Death (QBD) process. Using traditional results related to QBD processes, equilibrium performance measures are determined. More interestingly, we demonstrate the tractability of the computation of transient probability distribution of a QBD process. In the research literature, transient analysis of the QBD process was shown to be tractable in the Laplace transform domain. But in this research paper, we prove that the matrix exponential eQt arising in transient analysis (where Q is the generator matrix of the QBD process) can be computed directly in the time domain rendering efficient transient analysis of QBD process. Using the transient probability mass function of queue length, estimation of transient performance measures such as expected queue length, average delay, and tail distribution can be determined. Further, optimal adaptive routing algorithms for coupled queues can be designed.

Α new mixed δ-shock model with a change in shock distribution.

In this paper, reliability properties of a system that is subject to a sequence of shocks are investigated under a particular new change point model. According to the model, a change in the distribution of the shock magnitudes occurs upon the occurrence of a shock that is above a certain critical level. The system fails when the time between successive shocks is less than a given threshold, or the magnitude of a single shock is above a critical threshold. The survival function of the system is studied under both cases when the times between shocks follow discrete distribution and when the times between shocks follow continuous distribution. Matrix-based expressions are obtained for matrix-geometric discrete intershock times and for matrix-exponential continuous intershock times, as well.

An M/PH/1 queue with workload-dependent processing speed and vacations.

Motivated by the trade-off issue between delay performance and energy consumption in modern computer and communication systems, we consider a single-server queue with phase-type service requirements and with the following two special features: Firstly, the service speed is a piecewise constant function of the workload. Secondly, the server switches off when the system becomes empty, only to be activated again when the workload reaches a certain threshold. For this system, we obtain the steady-state workload distribution and its moments of any order. We use this result to choose the activation threshold such that a certain cost function, involving processing costs, activation costs and mean workload, is minimized.

Spatial inequalities of COVID-19 mortality rate in relation to socioeconomic and environmental factors across England.

In this study, we aimed to examine spatial inequalities of COVID-19 mortality rate in relation to spatial inequalities of socioeconomic and environmental factors across England. Specifically, we first explored spatial patterns of COVID-19 mortality rate in comparison to non-COVID-19 mortality rate. Subsequently, we established models to investigate contributions of socioeconomic and environmental factors to spatial variations of COVID-19 mortality rate across England (N = 317). Two newly developed specifications of spatial regression models were established successfully to estimate COVID-19 mortality rate (R2 = 0.49 and R2 = 0.793). The level of spatial inequalities of COVID-19 mortality is higher than that of non-COVID-19 mortality in England. Although global spatial association of COVID-19 mortality and non-COVID-19 mortality is positive, local spatial association of COVID-19 mortality and non-COVID-19 mortality is negative in some areas. Expectedly, hospital accessibility is negatively related to COVID-19 mortality rate. Percent of Asians, percent of Blacks, and unemployment rate are positively related to COVID-19 mortality rate. More importantly, relative humidity is negatively related to COVID-19 mortality rate. Moreover, among the spatial models estimated, the 'random effects specification of eigenvector spatial filtering model' outperforms the 'matrix exponential spatial specification of spatial autoregressive model'.

Matrix exponential based discriminant locality preserving projections for feature extraction.

Discriminant locality preserving projections (DLPP), which has shown good performances in pattern recognition, is a feature extraction algorithm based on manifold learning. However, DLPP suffers from the well-known small sample size (SSS) problem, where the number of samples is less than the dimension of samples. In this paper, we propose a novel matrix exponential based discriminant locality preserving projections (MEDLPP). The proposed MEDLPP method can address the SSS problem elegantly since the matrix exponential of a symmetric matrix is always positive definite. Nevertheless, the computational complexity of MEDLPP is high since it needs to solve a large matrix exponential eigenproblem. Then, in this paper, we also present an efficient algorithm for solving MEDLPP. Besides, the main idea for solving MEDLPP efficiently is also generalized to other matrix exponential based methods. The experimental results on some data sets demonstrate the proposed algorithm outperforms many state-of-the-art discriminant analysis methods.

Is a matrix exponential specification suitable for the modeling of spatial correlation structures?

This paper investigates the adequacy of the matrix exponential spatial specifications (MESS) as an alternative to the widely used spatial autoregressive models (SAR). To provide as complete a picture as possible, we extend the analysis to all the main spatial models governed by matrix exponentials comparing them with their spatial autoregressive counterparts. We propose a new implementation of Bayesian parameter estimation for the MESS model with vague prior distributions, which is shown to be precise and computationally efficient. Our implementations also account for spatially lagged regressors. We further allow for location-specific heterogeneity, which we model by including spatial splines. We conclude by comparing the performances of the different model specifications in applications to a real data set and by running simulations. Both the applications and the simulations suggest that the spatial splines are a flexible and efficient way to account for spatial heterogeneities governed by unknown mechanisms.

A computationally efficient algorithm for fitting ion channel parameters.

Continuous time Markov models have been widely used to describe ion channel kinetics, providing explicit representation of channel states and transitions. Fitting models to experimental data remains a computationally demanding task largely due to the high cost of model evaluation. Here, we propose a method to efficiently optimize model parameters and structure. Voltage clamp channel protocols can be decomposed into a series of fixed steps of constant voltage resulting in a set of linear systems of differential equations. Given the linear systems, ODE integration can be swapped for the faster matrix exponential routine. With our parallelized implementation, optimized models are able to reproduce a wide range of experimentally collected data within one minute, a 50 times speedup over ODE integration. •The cost of the objective function is reduced by employing the matrix exponential•The likelihood of convergence is improved by applying synchronous start simulated annealing•The approach was tested by optimizing parameters for a model of the cardiac voltage-gated Na+ channel, NaV1.5, and the KCNQ1 K+ channel.

Series: Utilization of Differential Equations and Methods for Solving Them in Medical Physics (3).

In this issue, simultaneous differential equations were introduced. These differential equations are often used in the field of medical physics. The methods for solving them were also introduced, which include Laplace transform and matrix methods. Some examples were also introduced, in which Laplace transform and matrix methods were applied to solving simultaneous differential equations derived from a three-compartment kinetic model for analyzing the glucose metabolism in tissues and Bloch equations for describing the behavior of the macroscopic magnetization in magnetic resonance imaging.In the next (final) issue, partial differential equations and various methods for solving them will be introduced together with some examples in medical physics.

Solving the chemical master equation by a fast adaptive finite state projection based on the stochastic simulation algorithm.

The mathematical framework of the chemical master equation (CME) uses a Markov chain to model the biochemical reactions that are taking place within a biological cell. Computing the transient probability distribution of this Markov chain allows us to track the composition of molecules inside the cell over time, with important practical applications in a number of areas such as molecular biology or medicine. However the CME is typically difficult to solve, since the state space involved can be very large or even countably infinite. We present a novel way of using the stochastic simulation algorithm (SSA) to reduce the size of the finite state projection (FSP) method. Numerical experiments that demonstrate the effectiveness of the reduction are included.