Fides: Reliable trust-region optimization for parameter estimation of ordinary differential equation models

Ordinary differential equation (ODE) models are widely used to study biochemical reactions in cellular networks since they effectively describe the temporal evolution of these networks using mass action kinetics. The parameters of these models are rarely known a priori and must instead be estimated by calibration using experimental data. Optimization-based calibration of ODE models on is often challenging, even for low-dimensional problems. Multiple hypotheses have been advanced to explain why biochemical model calibration is challenging, including non-identifiability of model parameters, but there are few comprehensive studies that test these hypotheses, likely because tools for performing such studies are also lacking. Nonetheless, reliable model calibration is essential for uncertainty analysis, model comparison, and biological interpretation.We implemented an established trust-region method as a modular Python framework (fides) to enable systematic comparison of different approaches to ODE model calibration involving a variety of Hessian approximation schemes. We evaluated fides on a recently developed corpus of biologically realistic benchmark problems for which real experimental data are available. Unexpectedly, we observed high variability in optimizer performance among different implementations of the same mathematical instructions (algorithms). Analysis of possible sources of poor optimizer performance identified limitations in the widely used Gauss-Newton, BFGS and SR1 Hessian approximation schemes. We addressed these drawbacks with a novel hybrid Hessian approximation scheme that enhances optimizer performance and outperforms existing hybrid approaches. When applied to the corpus of test models, we found that fides was on average more reliable and efficient than existing methods using a variety of criteria. We expect fides to be broadly useful for ODE constrained optimization problems in biochemical models and to be a foundation for future methods development.     

On the inference of complex phylogenetic networks by Markov Chain Monte-Carlo

For various species, high quality sequences and complete genomes are nowadays available for many individuals. This makes data analysis challenging, as methods need not only to be accurate, but also time efficient given the tremendous amount of data to process. In this article, we introduce an efficient method to infer the evolutionary history of individuals under the multispecies coalescent model in networks (MSNC). Phylogenetic networks are an extension of phylogenetic trees that can contain reticulate nodes, which allow to model complex biological events such as horizontal gene transfer, hybridization and introgression. We present a novel way to compute the likelihood of biallelic markers sampled along genomes whose evolution involved such events. This likelihood computation is at the heart of a Bayesian network inference method called SnappNet, as it extends the Snapp method inferring evolutionary trees under the multispecies coalescent model, to networks. SnappNet is available as a package of the well-known beast 2 software.Recently, the MCMC_BiMarkers method, implemented in PhyloNet, also extended Snapp to networks. Both methods take biallelic markers as input, rely on the same model of evolution and sample networks in a Bayesian framework, though using different methods for computing priors. However, SnappNet relies on algorithms that are exponentially more time-efficient on non-trivial networks. Using simulations, we compare performances of SnappNet and MCMC_BiMarkers. We show that both methods enjoy similar abilities to recover simple networks, but SnappNet is more accurate than MCMC_BiMarkers on more complex network scenarios. Also, on complex networks, SnappNet is found to be extremely faster than MCMC_BiMarkers in terms of time required for the likelihood computation. We finally illustrate SnappNet performances on a rice data set. SnappNet infers a scenario that is consistent with previous results and provides additional understanding of rice evolution.     

GPRuler: Metabolic gene-protein-reaction rules automatic reconstruction

Metabolic network models are increasingly being used in health care and industry. As a consequence, many tools have been released to automate their reconstruction process de novo. In order to enable gene deletion simulations and integration of gene expression data, these networks must include gene-protein-reaction (GPR) rules, which describe with a Boolean logic relationships between the gene products (e.g., enzyme isoforms or subunits) associated with the catalysis of a given reaction. Nevertheless, the reconstruction of GPRs still remains a largely manual and time consuming process. Aiming at fully automating the reconstruction process of GPRs for any organism, we propose the open-source python-based framework GPRuler. By mining text and data from 9 different biological databases, GPRuler can reconstruct GPRs starting either from just the name of the target organism or from an existing metabolic model. The performance of the developed tool is evaluated at small-scale level for a manually curated metabolic model, and at genome-scale level for three metabolic models related to Homo sapiens and Saccharomyces cerevisiae organisms. By exploiting these models as benchmarks, the proposed tool shown its ability to reproduce the original GPR rules with a high level of accuracy. In all the tested scenarios, after a manual investigation of the mismatches between the rules proposed by GPRuler and the original ones, the proposed approach revealed to be in many cases more accurate than the original models. By complementing existing tools for metabolic network reconstruction with the possibility to reconstruct GPRs quickly and with a few resources, GPRuler paves the way to the study of context-specific metabolic networks, representing the active portion of the complete network in given conditions, for organisms of industrial or biomedical interest that have not been characterized metabolically yet.     

A scale-free structure prior for graphical models with applications in functional genomics

The problem of reconstructing large-scale, gene regulatory networks from gene expression data has garnered considerable attention in bioinformatics over the past decade with the graphical modeling paradigm having emerged as a popular framework for inference. Analysis in a full Bayesian setting is contingent upon the assignment of a so-called structure prior-a probability distribution on networks, encoding a priori biological knowledge either in the form of supplemental data or high-level topological features. A key topological consideration is that a wide range of cellular networks are approximately scale-free, meaning that the fraction, , of nodes in a network with degree is roughly described by a power-law with exponent between and . The standard practice, however, is to utilize a random structure prior, which favors networks with binomially distributed degree distributions. In this paper, we introduce a scale-free structure prior for graphical models based on the formula for the probability of a network under a simple scale-free network model. Unlike the random structure prior, its scale-free counterpart requires a node labeling as a parameter. In order to use this prior for large-scale network inference, we design a novel Metropolis-Hastings sampler for graphical models that includes a node labeling as a state space variable. In a simulation study, we demonstrate that the scale-free structure prior outperforms the random structure prior at recovering scale-free networks while at the same time retains the ability to recover random networks. We then estimate a gene association network from gene expression data taken from a breast cancer tumor study, showing that scale-free structure prior recovers hubs, including the previously unknown hub SLC39A6, which is a zinc transporter that has been implicated with the spread of breast cancer to the lymph nodes. Our analysis of the breast cancer expression data underscores the value of the scale-free structure prior as an instrument to aid in the identification of candidate hub genes with the potential to direct the hypotheses of molecular biologists, and thus drive future experiments.     

A condition for cooperation in a game on complex networks

We study a condition of favoring cooperation in Prisoner's Dilemma game on complex networks. There are two kinds of players: cooperators and defectors. Cooperators pay a benefit b to their neighbors at a cost c, whereas defectors only receive a benefit. The game is a death-birth process with weak selection. Although it has been widely thought that b/c>〈k〉 is a condition of favoring cooperation (Ohtsuki et al., 2006), we find that b/c>〈k_nn〉 is the condition. We also show that among three representative networks, namely, regular, random, and scale-free, a regular network favors cooperation the most, whereas a scale-free network favors cooperation the least. In an ideal scale-free network, cooperation is never realized. Whether or not the scale-free network and network heterogeneity favor cooperation depends on the details of the game, although it is occasionally believed that these favor cooperation irrespective of the game structure.     

A mathematical model for the distribution of hemodynamic parameters in the human retinal microvascular network

To quantitatively assess the arteriovenous distribution of hemodynamic parameters throughout the microvascular network of the human retina, we constructed a retinal microcirculatory model consisting of a dichotomous symmetric branching system. This system is characterized by a diameter exponent of 2.85, instead of 3 as dictated by Murray’s law, except for the capillary networks. The value of 2.85 was the sum of a fractal dimension (1.70) and a branch exponent (1.15) of the retinal vasculature. Following the feeding artery (central retinal artery), each bifurcation was recursively developed at a distance of an individual branch length [L(r) = 7.4r ^1.15] by a centrifugal scheme. The venular tree was formed in the same way. Using this model, we evaluated hemodynamic parameters, including blood pressure, blood flow, blood velocity, shear rate, and shear stress, within the retinal microcirculatory network as a function of vessel diameter. The arteriovenous distributions of blood pressure and velocity in the simulation were consistent with in vivo measurements in the human retina and other vascular beds of small animals. We therefore conclude that the current theoretical model was useful for quantifying hemodynamics as a function of vessel diameter within the retinal microvascular network.     

Fast Generation of Sparse Random Kernel Graphs

The development of kernel-based inhomogeneous random graphs has provided models that are flexible enough to capture many observed characteristics of real networks, and that are also mathematically tractable. We specify a class of inhomogeneous random graph models, called random kernel graphs, that produces sparse graphs with tunable graph properties, and we develop an efficient generation algorithm to sample random instances from this model. As real-world networks are usually large, it is essential that the run-time of generation algorithms scales better than quadratically in the number of vertices n. We show that for many practical kernels our algorithm runs in time at most 𝒪(n(logn)²). As a practical example we show how to generate samples of power-law degree distribution graphs with tunable assortativity.     

Intrinsic properties of Boolean dynamics in complex networks

We study intrinsic properties of attractor in Boolean dynamics of complex networks with scale-free topology, comparing with those of the so-called Kauffman's random Boolean networks. We numerically study both frozen and relevant nodes in each attractor in the dynamics of relatively small networks (20?N?200). We investigate numerically robustness of an attractor to a perturbation. An attractor with cycle length of ?_c in a network of size N consists of ?_c states in the state space of 2^N states; each attractor has the arrangement of N nodes, where the cycle of attractor sweeps ?_c states. We define a perturbation as a flip of the state on a single node in the attractor state at a given time step. We show that the rate between unfrozen and relevant nodes in the dynamics of a complex network with scale-free topology is larger than that in Kauffman's random Boolean network model. Furthermore, we find that in a complex scale-free network with fluctuation of the in-degree number, attractors are more sensitive to a state flip for a highly connected node (i.e. input-hub node) than to that for a less connected node. By some numerical examples, we show that the number of relevant nodes increases, when an input-hub node is coincident with and/or connected with an output-hub node (i.e. a node with large output-degree) one another.     

Emergence of Super Cooperation of Prisoner’s Dilemma Games on Scale-Free Networks

Recently, the authors proposed a quantum prisoner’s dilemma game based on the spatial game of Nowak and May, and showed that the game can be played classically. By using this idea, we proposed three generalized prisoner’s dilemma (GPD, for short) games based on the weak Prisoner’s dilemma game, the full prisoner’s dilemma game and the normalized Prisoner’s dilemma game, written by GPD_W, GPD_F and GPD_N respectively. Our games consist of two players, each of which has three strategies: cooperator (C), defector (D) and super cooperator (denoted by Q), and have a parameter γ to measure the entangled relationship between the two players. We found that our generalised prisoner’s dilemma games have new Nash equilibrium principles, that entanglement is the principle of emergence and convergence (i.e., guaranteed emergence) of super cooperation in evolutions of our generalised prisoner’s dilemma games on scale-free networks, that entanglement provides a threshold for a phase transition of super cooperation in evolutions of our generalised prisoner’s dilemma games on scale-free networks, that the role of heterogeneity of the scale-free networks in cooperations and super cooperations is very limited, and that well-defined structures of scale-free networks allow coexistence of cooperators and super cooperators in the evolutions of the weak version of our generalised prisoner’s dilemma games.     

Asymptotic states and topological structure of an activation-deactivation chemical network

The influence of the topology on the asymptotic states of a network of interacting chemical species has been studied by simulating its time evolution. Random and scale-free networks have been designed to support relevant features of activation-deactivation reactions networks (mapping signal transduction networks) and the system of ordinary differential equations associated to the dynamics has been numerically solved. We analysed stationary states of the dynamics as a function of the network's connectivity and of the distribution of the chemical species on the network; we found important differences between the two topologies in the regime of low connectivity. In particular, only for low connected scale-free networks it is possible to find zero activity patterns as stationary states of the dynamics which work as signal off-states. Asymptotic features of random and scale-free networks become similar as the connectivity increases.     

Network-State Modulation of Power-Law Frequency-Scaling in Visual Cortical Neurons

Various types of neural-based signals, such as EEG, local field potentials and intracellular synaptic potentials, integrate multiple sources of activity distributed across large assemblies. They have in common a power-law frequency-scaling structure at high frequencies, but it is still unclear whether this scaling property is dominated by intrinsic neuronal properties or by network activity. The latter case is particularly interesting because if frequency-scaling reflects the network state it could be used to characterize the functional impact of the connectivity. In intracellularly recorded neurons of cat primary visual cortex in vivo, the power spectral density of V_m activity displays a power-law structure at high frequencies with a fractional scaling exponent. We show that this exponent is not constant, but depends on the visual statistics used to drive the network. To investigate the determinants of this frequency-scaling, we considered a generic recurrent model of cortex receiving a retinotopically organized external input. Similarly to the in vivo case, our in computo simulations show that the scaling exponent reflects the correlation level imposed in the input. This systematic dependence was also replicated at the single cell level, by controlling independently, in a parametric way, the strength and the temporal decay of the pairwise correlation between presynaptic inputs. This last model was implemented in vitro by imposing the correlation control in artificial presynaptic spike trains through dynamic-clamp techniques. These in vitro manipulations induced a modulation of the scaling exponent, similar to that observed in vivo and predicted in computo. We conclude that the frequency-scaling exponent of the V_m reflects stimulus-driven correlations in the cortical network activity. Therefore, we propose that the scaling exponent could be used to read-out the “effective” connectivity responsible for the dynamical signature of the population signals measured at different integration levels, from Vm to LFP, EEG and fMRI.     

Characterization of the Complete Mitochondrial Genome of Cerura menciana and Comparison with Other Lepidopteran Insects

The complete mitochondrial genome (mitogenome) of Cerura menciana (Lepidoptera: Notodontidae) was sequenced and analyzed in this study. The mitogenome is a circular molecule of 15,369 bp, containing 13 protein-coding genes (PCGs), two ribosomal RNA (rRNA) genes, 22 transfer RNA (tRNA) genes and a A+T-rich region. The positive AT skew (0.031) indicated that more As than Ts were present. All PCGs were initiated by ATN codons, except for the cytochrome c oxidase subunit 1 (cox1) gene, which was initiated by CAG. Two of the 13 PCGs contained the incomplete termination codon T or TA, while the others were terminated with the stop codon TAA. The A+T-rich region was 372 bp in length and consisted of an ‘ATAGA’ motif followed by an 18 bp poly-T stretch, a microsatellite-like (AT)₈ and a poly-A element upstream of the trnM gene. Results examining codon usage indicated that Asn, Ile, Leu2, Lys, Tyr and Phe were the six most frequently occurring amino acids, while Cys was the rarest. Phylogenetic relationships, analyzed based on the nucleotide sequences of the 13 PCGs from other insect mitogenomes, confirmed that C. menciana belongs to the Notodontidae family.     

Spatial network structure and metapopulation persistence

We explore the relationship between network structure and dynamics by relating the topology of spatial networks with its underlying metapopulation abundance. Metapopulation abundance is largely affected by the architecture of the spatial network, although this effect depends on demographic parameters here represented by the extinction-to-colonization ratio (e/c). Thus, for moderate to large e/c-values, regional abundance grows with the heterogeneity of the network, with uniform or random networks having the lowest regional abundances, and scale-free networks having the largest abundance. However, the ranking is reversed for low extinction probabilities, with heterogeneous networks showing the lowest relative abundance. We further explore the mechanisms underlying such results by relating a node's incidence (average number of time steps the node is occupied) with its degree, and with the average degree of the nodes it interacts with. These results demonstrate the importance of spatial network structure to understanding metapopulation abundance, and serve to determine under what circumstances information on network structure should be complemented with information on the species life-history traits to understand persistence in heterogeneous environments.     

Ordering Dynamics in Neuron Activity Pattern Model: An Insight to Brain Functionality

We study the domain ordering kinetics in d = 2 ferromagnets which corresponds to populated neuron activities with both long-ranged interactions, V(r) ∼ r ⁻ⁿ and short-ranged interactions. We present the results from comprehensive Monte Carlo (MC) simulations for the nonconserved Ising model with n ≥ 2, interaction range considering near and far neighbors. Our model results could represent the long-ranged neuron kinetics (n ≤ 4) in consistent with the same dynamical behaviour of short-ranged case (n ≥ 4) at far below and near criticality. We found that emergence of fast and slow kinetics of long and short ranged case could imitate the formation of connections among near and distant neurons. The calculated characteristic length scale in long-ranged interaction is found to be n independent (L(t) ∼ t ^1/(n−2)), whereas short-ranged interaction follows L(t) ∼ t ^1/2 law and approximately preserve universality in domain kinetics. Further, we did the comparative study of phase ordering near the critical temperature which follows different behaviours of domain ordering near and far critical temperature but follows universal scaling law.     

Scaling and structure of dicotyledonous leaf venation networks

There have been numerous attempts to derive general models for the structure and function of resource delivery networks in biology. Such theories typically predict the quantitative structure of vascular networks across scales. For example, fractal branching models of plant structure predict that the network dimensions within plant stems or leaves should be scale-free. However, very few empirical examples of such networks are available with which to evaluate such hypotheses. Here, we apply recently developed leaf network extraction software to a global leaf dataset. We find that leaf networks are neither entirely scale-free nor governed entirely by a characteristic scale. Indeed, we find many network properties, such as vein length distributions, which are governed by characteristic scales, and other network properties, notably vein diameter distributions, which are typified by power-law behaviour. Our findings suggest that theories of network structure will remain incomplete until they address the multiple constraints on network architecture.     

Identification procedure in a model of single fibre action potential – Part II: Global approach and experimental results

The present paper describes a global procedure for estimating all the synthesis parameters that generate a single fibre action potential (SFAP) in the Dimitrov–Dimitrova (D–D) convolutional model. We call this inverse problem Identification Procedure, and it is presented in two parts, this paper being the second. The procedure incorporates the candidate pair (CP) method developed in Part I, which provides the values of radial distance r and fibre diameter d of the simulated SFAP that best matches a potential under study. The CP-method required prior knowledge of all the excitation parameters. However, since the Identification Procedure makes no assumption about the excitation, multiple combinations of the synthesis parameters result in very similar SFAPs whose shape is close the signal under study. Analysis of the possible combinations reveals that r and d can be modelled as two jointly Gaussian random variables. The interest of the Identification Procedure is that, for a certain SFAP, it provides estimates of r and d, along with estimates of different parameters that determine the IAP waveform. Moreover, the procedure is able to determine the degree of error that accompanies the estimation of r and d.