Meta-models as a straightforward approach to the sensitivity analysis of complex models

Complex simulation models are important tools in applied ecological and conservation research. However sensitivity analysis of this important class of models can be difficult to conduct. High level interactions and non-linear responses are common in complex simulations, and this necessitates a global sensitivity analysis, where each parameter is tested at a range of values, and in combination with changes in many other parameters. We reviewed the literature, searching for population viability analyses that used simulation models. We found only 9 out of the 122 simulation population viability analysis used global sensitivity analysis. This result is typical of other simulation models in applied ecology, where global sensitivity analysis is rare. We then demonstrate how to conduct a meta-modeling sensitivity analysis, where a simpler statistically fit function (the meta-model, also known as the surrogate model or emulator) is used to approximate the behavior of the complicated simulation. This simpler meta-model is interrogated to inform on the behavior of simulation model. We fit two example meta-models, a generalized linear model and a boosted regression tree, to exemplify the approach. Our hope is that by going through these techniques thoroughly they will become more widely adopted.     

Rapid Global Fitting of Large Fluorescence Lifetime Imaging Microscopy Datasets

Fluorescence lifetime imaging (FLIM) is widely applied to obtain quantitative information from fluorescence signals, particularly using Förster Resonant Energy Transfer (FRET) measurements to map, for example, protein-protein interactions. Extracting FRET efficiencies or population fractions typically entails fitting data to complex fluorescence decay models but such experiments are frequently photon constrained, particularly for live cell or in vivo imaging, and this leads to unacceptable errors when analysing data on a pixel-wise basis. Lifetimes and population fractions may, however, be more robustly extracted using global analysis to simultaneously fit the fluorescence decay data of all pixels in an image or dataset to a multi-exponential model under the assumption that the lifetime components are invariant across the image (dataset). This approach is often considered to be prohibitively slow and/or computationally expensive but we present here a computationally efficient global analysis algorithm for the analysis of time-correlated single photon counting (TCSPC) or time-gated FLIM data based on variable projection. It makes efficient use of both computer processor and memory resources, requiring less than a minute to analyse time series and multiwell plate datasets with hundreds of FLIM images on standard personal computers. This lifetime analysis takes account of repetitive excitation, including fluorescence photons excited by earlier pulses contributing to the fit, and is able to accommodate time-varying backgrounds and instrument response functions. We demonstrate that this global approach allows us to readily fit time-resolved fluorescence data to complex models including a four-exponential model of a FRET system, for which the FRET efficiencies of the two species of a bi-exponential donor are linked, and polarisation-resolved lifetime data, where a fluorescence intensity and bi-exponential anisotropy decay model is applied to the analysis of live cell homo-FRET data. A software package implementing this algorithm, FLIMfit, is available under an open source licence through the Open Microscopy Environment.     

Predictions of response to temperature are contingent on model choice and data quality

The equations used to account for the temperature dependence of biological processes, including growth and metabolic rates, are the foundations of our predictions of how global biogeochemistry and biogeography change in response to global climate change. We review and test the use of 12 equations used to model the temperature dependence of biological processes across the full range of their temperature response, including supra‐ and suboptimal temperatures. We focus on fitting these equations to thermal response curves for phytoplankton growth but also tested the equations on a variety of traits across a wide diversity of organisms. We found that many of the surveyed equations have comparable abilities to fit data and equally high requirements for data quality (number of test temperatures and range of response captured) but lead to different estimates of cardinal temperatures and of the biological rates at these temperatures. When these rate estimates are used for biogeographic predictions, differences between the estimates of even the best‐fitting models can exceed the global biological change predicted for a decade of global warming. As a result, studies of the biological response to global changes in temperature must make careful consideration of model selection and of the quality of the data used for parametrizing these models.     

Non-linear analysis of GeneChip arrays

The application of microarray hybridization theory to Affymetrix GeneChip data has been a recent focus for data analysts. It has been shown that the hyperbolic Langmuir isotherm captures the shape of the signal response to concentration of Affymetrix GeneChips. We demonstrate that existing linear fit methods for extracting gene expression measures are not well adapted for the effect of saturation resulting from surface adsorption processes. In contrast to the most popular methods, we fit background and concentration parameters within a single global fitting routine instead of estimating the background before obtaining gene expression measures. We describe a non-linear multi-chip model of the perfect match signal that effectively allows for the separation of specific and non-specific components of the microarray signal and avoids saturation bias in the high-intensity range. Multimodel inference, incorporated within the fitting routine, allows a quantitative selection of the model that best describes the observed data. The performance of this method is evaluated on publicly available datasets, and comparisons to popular algorithms are presented.     

Propagation of local interactions create global gap structure and dynamics in a tropical rainforest

Gap dynamics in tropical forests are of interest because an understanding of them can help to predict canopy structure and biodiversity. We present a simple cellular automaton model that is capable of capturing many of the trends seen in the canopy gap pattern of a complex tropical rainforest on the Barro Colorado Island (BCI) using a single set of model parameters. We fit the global and local densities, the cluster size distributions, and two correlation functions, for gaps, gap formations, and gap closures determined from a spatial map of the forest (1983-1984). To the best of our knowledge, this is the first report that the cluster size distributions of gap formations and closures in the BCI are both power laws. An important element in the model is that when a transition from gap to non-gap (closure), or vice versa (formation), occurs, this transition is allowed to expand into adjacent cells in order to make different cluster sizes of transitions. Model results are in excellent agreement with reported field data. The propagation of local interactions is necessary in order to obtain the complex dynamics of the gap pattern. We also establish a connection between the global and local densities via the neighborhood-dependent transition rates and the effective global transition rates.     

Temperature Dependence of Soil Respiration Modulated by Thresholds in Soil Water Availability Across European Shrubland Ecosystems

Soil respiration (SR) is a major component of the global carbon cycle and plays a fundamental role in ecosystem feedback to climate change. Empirical modelling is an essential tool for predicting ecosystem responses to environmental change, and also provides important data for calibrating and corroborating process-based models. In this study, we evaluated the performance of three empirical temperature–SR response functions (exponential, Lloyd–Taylor and Gaussian) at seven shrublands located within three climatic regions (Atlantic, Mediterranean and Continental) across Europe. We investigated the performance of SR models by including the interaction between soil moisture and soil temperature. We found that the best fit for the temperature functions depended on the site-specific climatic conditions. Including soil moisture, we identified thresholds in the three different response functions that improved the model fit in all cases. The direct soil moisture effect on SR, however, was weak at the annual time scale. We conclude that the exponential soil temperature function may only be a good predictor for SR in a narrow temperature range, and that extrapolating predictions for future climate based on this function should be treated with caution as modelled outputs may underestimate SR. The addition of soil moisture thresholds improved the model fit at all sites, but had a far greater ecological significance in the wet Atlantic shrubland where a fundamental change in the soil CO₂ efflux would likely have an impact on the whole carbon budget.     

Multi-scale inference of interaction rules in animal groups using Bayesian model selection

Inference of interaction rules of animals moving in groups usually relies on an analysis of large scale system behaviour. Models are tuned through repeated simulation until they match the observed behaviour. More recent work has used the fine scale motions of animals to validate and fit the rules of interaction of animals in groups. Here, we use a Bayesian methodology to compare a variety of models to the collective motion of glass prawns (Paratya australiensis). We show that these exhibit a stereotypical 'phase transition', whereby an increase in density leads to the onset of collective motion in one direction. We fit models to this data, which range from: a mean-field model where all prawns interact globally; to a spatial Markovian model where prawns are self-propelled particles influenced only by the current positions and directions of their neighbours; up to non-Markovian models where prawns have 'memory' of previous interactions, integrating their experiences over time when deciding to change behaviour. We show that the mean-field model fits the large scale behaviour of the system, but does not capture fine scale rules of interaction, which are primarily mediated by physical contact. Conversely, the Markovian self-propelled particle model captures the fine scale rules of interaction but fails to reproduce global dynamics. The most sophisticated model, the non-Markovian model, provides a good match to the data at both the fine scale and in terms of reproducing global dynamics. We conclude that prawns' movements are influenced by not just the current direction of nearby conspecifics, but also those encountered in the recent past. Given the simplicity of prawns as a study system our research suggests that self-propelled particle models of collective motion should, if they are to be realistic at multiple biological scales, include memory of previous interactions and other non-Markovian effects.     

When the optimal is not the best: parameter estimation in complex biological models

Background

The vast computational resources that became available during the past decade enabled the development and simulation of increasingly complex mathematical models of cancer growth. These models typically involve many free parameters whose determination is a substantial obstacle to model development. Direct measurement of biochemical parameters in vivo is often difficult and sometimes impracticable, while fitting them under data-poor conditions may result in biologically implausible values.

Results

We discuss different methodological approaches to estimate parameters in complex biological models. We make use of the high computational power of the Blue Gene technology to perform an extensive study of the parameter space in a model of avascular tumor growth. We explicitly show that the landscape of the cost function used to optimize the model to the data has a very rugged surface in parameter space. This cost function has many local minima with unrealistic solutions, including the global minimum corresponding to the best fit.

Conclusions

The case studied in this paper shows one example in which model parameters that optimally fit the data are not necessarily the best ones from a biological point of view. To avoid force-fitting a model to a dataset, we propose that the best model parameters should be found by choosing, among suboptimal parameters, those that match criteria other than the ones used to fit the model. We also conclude that the model, data and optimization approach form a new complex system and point to the need of a theory that addresses this problem more generally.     

SANS/SAXS study of the BSA solvation properties in aqueous urea solutions via a global fit approach

We report on the solvation properties and intermolecular interactions of a model protein (bovine serum albumine, BSA) in urea aqueous solutions, as obtained by combining small-angle neutron and X-ray scattering experiments. According to a global fit strategy, all the whole set of scattering curves are analysed by considering a unique model which includes the BSA structure, the protein-protein interactions and the thermodynamic exchange process of water/urea molecules at the protein solvent interface. As a main result, the equilibrium constant that accounts for the difference in composition between the bulk solvent and the protein solvation layer is derived. Results confirm that urea preferentially sticks to the protein surface, inducing a noticeable change in both the repulsive and the attractive interaction potentials.     

Partitioning global change: Assessing the relative importance of changes in climate and land cover for changes in avian distribution

Understanding the relative impact of climate change and land cover change on changes in avian distribution has implications for the future course of avian distributions and appropriate management strategies. Due to the dynamic nature of climate change, our goal was to investigate the processes that shape species distributions, rather than the current distributional patterns. To this end, we analyzed changes in the distribution of Eastern Wood Pewees (Contopus virens) and Red‐eyed Vireos (Vireo olivaceus) from 1997 to 2012 using Breeding Bird Survey data and dynamic correlated‐detection occupancy models. We estimated the local colonization and extinction rates of these species in relation to changes in climate (hours of extreme temperature) and changes in land cover (amount of nesting habitat). We fit six nested models to partition the deviance explained by spatial and temporal components of land cover and climate. We isolated the temporal components of environmental variables because this is the essence of global change. For both species, model fit was significantly improved when we modeled vital rates as a function of spatial variation in climate and land cover. Model fit improved only marginally when we added temporal variation in climate and land cover to the model. Temporal variation in climate explained more deviance than temporal variation in land cover, although both combined only explained 20% (Eastern Wood Pewee) and 6% (Red‐eyed Vireo) of temporal variation in vital rates. Our results showing a significant correlation between initial occupancy and environmental covariates are consistent with biological expectation and previous studies. The weak correlation between vital rates and temporal changes in covariates indicated that we have yet to identify the most relevant components of global change influencing the distributions of these species and, more importantly, that spatially significant covariates are not necessarily driving temporal shifts in avian distributions.     

Optimal estimation of ion-channel kinetics from macroscopic currents

Markov modeling provides an effective approach for modeling ion channel kinetics. There are several search algorithms for global fitting of macroscopic or single-channel currents across different experimental conditions. Here we present a particle swarm optimization(PSO)-based approach which, when used in combination with golden section search (GSS), can fit macroscopic voltage responses with a high degree of accuracy (errors within 1%) and reasonable amount of calculation time (less than 10 hours for 20 free parameters) on a desktop computer. We also describe a method for initial value estimation of the model parameters, which appears to favor identification of global optimum and can further reduce the computational cost. The PSO-GSS algorithm is applicable for kinetic models of arbitrary topology and size and compatible with common stimulation protocols, which provides a convenient approach for establishing kinetic models at the macroscopic level.     

The binding of W7, an inhibitor of striated muscle contraction, to cardiac troponin C

W7 is a well-characterized calmodulin antagonist. It decreases the maximal tension and rate of ATP hydrolysis in cardiac muscle fibers. Cardiac troponin C (cTnC) has been previously implicated as the mechanistically significant target for W7 in the myofilament. Two-dimensional NMR spectra ({1H,15N}- and {1H,13C}-HSQCs) were used to monitor the Ca2+-dependent binding of W7 to cTnC. Titration of cTnC x 3Ca2+ with W7 indicated binding to both domains of the protein. We examined the binding of W7 to the separated domains of cTnC to simplify the spectral analysis. In the titration of the C-terminal domain (cCTnC x 2Ca2+), the spectral peaks originating from a subset of residues changed nonuniformly, and could not be well-described as single-site binding. A global fit of the cCTnC x 2Ca2+ titration data to a two-site, sequential binding model (47 residues simultaneously fit) yielded a dissociation constant (Kd1) of 0.85-0.91 mM for the singly bound state, with the second dissociation constant fit to 3.40-3.65 mM (> or = 4 x Kd1). The titration data for the N-terminal domain (cNTnC x Ca2+) was globally fit to single-site binding model with a Kd of 0.15-0.30 mM (41 residues fit). The data are consistent with W7 binding to each domain's major hydrophobic pocket, coordinating side chains responsible for liganding cTnI. When in muscle fibers, W7 may compete with cTnI for target sites on cTnC.     

The kinetics of affinity-mediated cell-surface attachment

Data and a semi-empirical model are presented that describe the affinity interaction of yeast cells with a Concanavalin A derivatised surface. The model uses 3 parameters to describe the time course of cell attachment from a flowing suspension of yeast cells, over a range of flow rates, and gives an effective global fit to the data obtained. Further modifications allow the effects of a soluble competitor (glucose) on binding to be quantified in terms of a saturation effect, and an effective global fit is obtained. A comparison was made between the relationship between steady-state attached fraction and applied shear with similar data reported earlier (Ming, F. et al, 1998) for the detachment of pre-adsorbed cells. This shows that there is an order of magnitude difference between the forces required to effect complete detachment in the two systems, and that the nature of the relationship between shear and attached fraction is profoundly different. The magnitude of this time-dependent stabilization might be explained in terms of a progressive reorientation of cell relative to the surface such that the number of bonds is maximized.     

Branching,blending, and the evolution of cultural similarities and differences among human populations

It has been claimed that blending processes such as trade and exchange have always been more important in the evolution of cultural similarities and differences among human populations than the branching process of population fissioning. In this paper, we report the results of a novel comparative study designed to shed light on this claim. We fitted the bifurcating tree model that biologists use to represent the relationships of species to 21 biological data sets that have been used to reconstruct the relationships of species and/or higher level taxa and to 21 cultural data sets. We then compared the average fit between the biological data sets and the model with the average fit between the cultural data sets and the model. Given that the biological data sets can be confidently assumed to have been structured by speciation, which is a branching process, our assumption was that, if cultural evolution is dominated by blending processes, the fit between the bifurcating tree model and the cultural data sets should be significantly worse than the fit between the bifurcating tree model and the biological data sets. Conversely, if cultural evolution is dominated by branching processes, the fit between the bifurcating tree model and the cultural data sets should be no worse than the fit between the bifurcating tree model and the biological data sets. We found that the average fit between the cultural data sets and the bifurcating tree model was not significantly different from the fit between the biological data sets and the bifurcating tree model. This indicates that the cultural data sets are not less tree-like than are the biological data sets. As such, our analysis does not support the suggestion that blending processes have always been more important than branching processes in cultural evolution. We conclude from this that, rather than deciding how cultural evolution has proceeded a priori, researchers need to ascertain which model or combination of models is relevant in a particular case and why.     

Mechanistic theory and modelling of complex food-web dynamics in Lake Constance

Mechanistic understanding of consumer-resource dynamics is critical to predicting the effects of global change on ecosystem structure, function and services. Such understanding is severely limited by mechanistic models' inability to reproduce the dynamics of multiple populations interacting in the field. We surpass this limitation here by extending general consumer-resource network theory to the complex dynamics of a specific ecosystem comprised by the seasonal biomass and production patterns in a pelagic food web of a large, well-studied lake. We parameterised our allometric trophic network model of 24 guilds and 107 feeding relationships using the lake's food web structure, initial spring biomasses and body-masses. Adding activity respiration, the detrital loop, minimal abiotic forcing, prey resistance and several empirically observed rates substantially increased the model's fit to the observed seasonal dynamics and the size-abundance distribution. This process illuminates a promising approach towards improving food-web theory and dynamic models of specific habitats.     

Global pattern and local variation in species–area relationships

Aim We conducted a meta‐analysis of species–area relationships (SARs) by combining several data sets and important covariates such as types of islands, taxonomic groups, latitude and spatial extent, in a hierarchical model framework to study global pattern and local variation in SARs and its consequences for prediction. Location One thousand nine hundred and eighteen islands from 94 SAR studies from around the world. Methods We developed a generalization of the power‐law SAR model, the HSARX model, which allows: (1) the inclusion of multiple focal parameters (intercept, slope, within‐study variance), (2) use of multiple effect modifiers based on a collection of SAR studies, and (3) modelling of the between‐ and within‐study variability. Results The global pattern in the SAR was the average of local SARs and had wide confidence intervals. The global SAR slope was 0.228 with 90% confidence limits of 0.059 and 0.412. The intercept, slope and within‐study variability of local SARs showed great heterogeneity as a result of the interaction of modifying covariates. Confidence intervals for these SAR parameters were narrower when other covariates in addition to area were accounted for, thus increasing the accuracy of the predictions for species richness. The significant effect of latitude and the interaction of latitude, taxa and island type on the SAR slope indicated that the ‘typical’ latitudinal diversity gradient can be reversed in isolated systems. Main conclusions The power‐law relationship underlying the HSARX model provides a good fit for non‐nested SARs across vastly different spatial scales by taking into account other covariates. The HSARX framework allows researchers to explore the complex interactions among SAR parameters and modifying variables, to explicitly study the scale dependence, and to make robust predictions on multiple levels (island, study, global) with associated prediction intervals. From a prediction perspective, it is not the global pattern but the local variation that matters.     

The gambin model provides a superior fit to species abundance distributions with a single free parameter: evidence,implementation and interpretation

The species abundance distribution (SAD) has been a central focus of community ecology for over fifty years, and is currently the subject of widespread renewed interest. The gambin model has recently been proposed as a model that provides a superior fit to commonly preferred SAD models. It has also been argued that the model's single parameter (α) presents a potentially informative ecological diversity metric, because it summarises the shape of the SAD in a single number. Despite this potential, few empirical tests of the model have been undertaken, perhaps because the necessary methods and software for fitting the model have not existed. Here, we derive a maximum likelihood method to fit the model, and use it to undertake a comprehensive comparative analysis of the fit of the gambin model. The functions and computational code to fit the model are incorporated in a newly developed free‐to‐download R package (gambin). We test the gambin model using a variety of datasets and compare the fit of the gambin model to fits obtained using the Poisson lognormal, logseries and zero‐sum multinomial distributions. We found that gambin almost universally provided a better fit to the data and that the fit was consistent for a variety of sample grain sizes. We demonstrate how α can be used to differentiate intelligibly between community structures of Azorean arthropods sampled in different land use types. We conclude that gambin presents a flexible model capable of fitting a wide variety of observed SAD data, while providing a useful index of SAD form in its single fitted parameter. As such, gambin has wide potential applicability in the study of SADs, and ecology more generally.     

A probabilistic methodology for integrating knowledge and experiments on biological networks.

Biological systems are traditionally studied by focusing on a specific subsystem, building an intuitive model for it, and refining the model using results from carefully designed experiments. Modern experimental techniques provide massive data on the global behavior of biological systems, and systematically using these large datasets for refining existing knowledge is a major challenge. Here we introduce an extended computational framework that combines formalization of existing qualitative models, probabilistic modeling, and integration of high-throughput experimental data. Using our methods, it is possible to interpret genomewide measurements in the context of prior knowledge on the system, to assign statistical meaning to the accuracy of such knowledge, and to learn refined models with improved fit to the experiments. Our model is represented as a probabilistic factor graph, and the framework accommodates partial measurements of diverse biological elements. We study the performance of several probabilistic inference algorithms and show that hidden model variables can be reliably inferred even in the presence of feedback loops and complex logic. We show how to refine prior knowledge on combinatorial regulatory relations using hypothesis testing and derive p-values for learned model features. We test our methodology and algorithms on a simulated model and on two real yeast models. In particular, we use our method to explore uncharacterized relations among regulators in the yeast response to hyper-osmotic shock and in the yeast lysine biosynthesis system. Our integrative approach to the analysis of biological regulation is demonstrated to synergistically combine qualitative and quantitative evidence into concrete biological predictions.