Hands-on parameter search for neural simulations by a MIDI-controller

Computational neuroscientists frequently encounter the challenge of parameter fitting--exploring a usually high dimensional variable space to find a parameter set that reproduces an experimental data set. One common approach is using automated search algorithms such as gradient descent or genetic algorithms. However, these approaches suffer several shortcomings related to their lack of understanding the underlying question, such as defining a suitable error function or getting stuck in local minima. Another widespread approach is manual parameter fitting using a keyboard or a mouse, evaluating different parameter sets following the users intuition. However, this process is often cumbersome and time-intensive. Here, we present a new method for manual parameter fitting. A MIDI controller provides input to the simulation software, where model parameters are then tuned according to the knob and slider positions on the device. The model is immediately updated on every parameter change, continuously plotting the latest results. Given reasonably short simulation times of less than one second, we find this method to be highly efficient in quickly determining good parameter sets. Our approach bears a close resemblance to tuning the sound of an analog synthesizer, giving the user a very good intuition of the problem at hand, such as immediate feedback if and how results are affected by specific parameter changes. In addition to be used in research, our approach should be an ideal teaching tool, allowing students to interactively explore complex models such as Hodgkin-Huxley or dynamical systems.     

A self-organizing state-space-model approach for parameter estimation in hodgkin-huxley-type models of single neurons

Traditional approaches to the problem of parameter estimation in biophysical models of neurons and neural networks usually adopt a global search algorithm (for example, an evolutionary algorithm), often in combination with a local search method (such as gradient descent) in order to minimize the value of a cost function, which measures the discrepancy between various features of the available experimental data and model output. In this study, we approach the problem of parameter estimation in conductance-based models of single neurons from a different perspective. By adopting a hidden-dynamical-systems formalism, we expressed parameter estimation as an inference problem in these systems, which can then be tackled using a range of well-established statistical inference methods. The particular method we used was Kitagawa's self-organizing state-space model, which was applied on a number of Hodgkin-Huxley-type models using simulated or actual electrophysiological data. We showed that the algorithm can be used to estimate a large number of parameters, including maximal conductances, reversal potentials, kinetics of ionic currents, measurement and intrinsic noise, based on low-dimensional experimental data and sufficiently informative priors in the form of pre-defined constraints imposed on model parameters. The algorithm remained operational even when very noisy experimental data were used. Importantly, by combining the self-organizing state-space model with an adaptive sampling algorithm akin to the Covariance Matrix Adaptation Evolution Strategy, we achieved a significant reduction in the variance of parameter estimates. The algorithm did not require the explicit formulation of a cost function and it was straightforward to apply on compartmental models and multiple data sets. Overall, the proposed methodology is particularly suitable for resolving high-dimensional inference problems based on noisy electrophysiological data and, therefore, a potentially useful tool in the construction of biophysical neuron models.     

Efficient parameter estimation for spatio-temporal models of pattern formation: case study of Drosophila melanogaster

MOTIVATION: Diffusable and non-diffusable gene products play a major role in body plan formation. A quantitative understanding of the spatio-temporal patterns formed in body plan formation, by using simulation models is an important addition to experimental observation. The inverse modelling approach consists of describing the body plan formation by a rule-based model, and fitting the model parameters to real observed data. In body plan formation, the data are usually obtained from fluorescent immunohistochemistry or in situ hybridizations. Inferring model parameters by comparing such data to those from simulation is a major computational bottleneck. An important aspect in this process is the choice of method used for parameter estimation. When no information on parameters is available, parameter estimation is mostly done by means of heuristic algorithms. RESULTS: We show that parameter estimation for pattern formation models can be efficiently performed using an evolution strategy (ES). As a case study we use a quantitative spatio-temporal model of the regulatory network for early development in Drosophila melanogaster. In order to estimate the parameters, the simulated results are compared to a time series of gene products involved in the network obtained with immunohistochemistry. We demonstrate that a (mu,lambda)-ES can be used to find good quality solutions in the parameter estimation. We also show that an ES with multiple populations is 5-140 times as fast as parallel simulated annealing for this case study, and that combining ES with a local search results in an efficient parameter estimation method.     

Universally sloppy parameter sensitivities in systems biology models

Quantitative computational models play an increasingly important role in modern biology. Such models typically involve many free parameters, and assigning their values is often a substantial obstacle to model development. Directly measuring in vivo biochemical parameters is difficult, and collectively fitting them to other experimental data often yields large parameter uncertainties. Nevertheless, in earlier work we showed in a growth-factor-signaling model that collective fitting could yield well-constrained predictions, even when it left individual parameters very poorly constrained. We also showed that the model had a “sloppy” spectrum of parameter sensitivities, with eigenvalues roughly evenly distributed over many decades. Here we use a collection of models from the literature to test whether such sloppy spectra are common in systems biology. Strikingly, we find that every model we examine has a sloppy spectrum of sensitivities. We also test several consequences of this sloppiness for building predictive models. In particular, sloppiness suggests that collective fits to even large amounts of ideal time-series data will often leave many parameters poorly constrained. Tests over our model collection are consistent with this suggestion. This difficulty with collective fits may seem to argue for direct parameter measurements, but sloppiness also implies that such measurements must be formidably precise and complete to usefully constrain many model predictions. We confirm this implication in our growth-factor-signaling model. Our results suggest that sloppy sensitivity spectra are universal in systems biology models. The prevalence of sloppiness highlights the power of collective fits and suggests that modelers should focus on predictions rather than on parameters.     

Fitting Membrane Resistance along with Action Potential Shape in Cardiac Myocytes Improves Convergence: Application of a Multi-Objective Parallel Genetic Algorithm

Fitting parameter sets of non-linear equations in cardiac single cell ionic models to reproduce experimental behavior is a time consuming process. The standard procedure is to adjust maximum channel conductances in ionic models to reproduce action potentials (APs) recorded in isolated cells. However, vastly different sets of parameters can produce similar APs. Furthermore, even with an excellent AP match in case of single cell, tissue behaviour may be very different. We hypothesize that this uncertainty can be reduced by additionally fitting membrane resistance (R_m). To investigate the importance of R_m, we developed a genetic algorithm approach which incorporated R_m data calculated at a few points in the cycle, in addition to AP morphology. Performance was compared to a genetic algorithm using only AP morphology data. The optimal parameter sets and goodness of fit as computed by the different methods were compared. First, we fit an ionic model to itself, starting from a random parameter set. Next, we fit the AP of one ionic model to that of another. Finally, we fit an ionic model to experimentally recorded rabbit action potentials. Adding the extra objective (R_m, at a few voltages) to the AP fit, lead to much better convergence. Typically, a smaller MSE (mean square error, defined as the average of the squared error between the target AP and AP that is to be fitted) was achieved in one fifth of the number of generations compared to using only AP data. Importantly, the variability in fit parameters was also greatly reduced, with many parameters showing an order of magnitude decrease in variability. Adding R_m to the objective function improves the robustness of fitting, better preserving tissue level behavior, and should be incorporated.     

A numerical analysis model for interpretation of flow cytometric studies of ex vivo phagocytosis

The study of ex vivo phagocytosis via flow cytometry requires that one distinguish experimentally between uptake and adsorption of fluorescently labeled targets by phagocytes. Removal of the latter quantity from the analysis is the most common means of analyzing such data. Because the probability of phagocytosis is a function of the probability of adsorption, and because partially quenched fluorescence after uptake often overlaps with that of negative controls, this approach is suboptimal at best. Here, we describe a numerical analysis model which overcomes these limitations. We posit that the random adsorption of targets to macrophages, and subsequent phagocytosis, is a function of three parameters: the ratio of targets to macrophages (m), the mean fluorescence intensity imparted to the phagocyte by the internalized target (alpha), and the probability of phagocytosis per adsorbed target (p). The potential values of these parameters define a parameter space and their values at any point in parameter space can be used to predict the fraction of adsorption(+) and [adsorption(-), phagocytosis(+)] cells that might be observed experimentally. By systematically evaluating the points in parameter space for the latter two values and comparing them to experimental data, the model arrives at sets of parameter values that optimally predict such data. Using activated THP-1 cells as macrophages and platelets as targets, we validate the model by demonstrating that it can distinguish between the effects of experimental changes in m, alpha, and p. Finally, we use the model to demonstrate that platelets from a congenitally thrombocytopenic WAS patient show an increased probability of ex vivo phagocytosis. This finding correlates with other evidence that rapid in vivo platelet consumption contributes significantly to the thrombocytopenia of WAS. Our numerical analysis method represents a useful and innovative approach to multivariate analysis.     

A simple work flow for biologically inspired model reduction - application to early JAK-STAT signaling

Background  

Modeling of biological pathways is a key issue in systems biology. When constructing a model, it is tempting to incorporate all known interactions of pathway species, which results in models with a large number of unknown parameters. Fortunately, unknown parameters need not necessarily be measured directly, but some parameter values can be estimated indirectly by fitting the model to experimental data. However, parameter fitting, or, more precisely, maximum likelihood parameter estimation, only provides valid results, if the complexity of the model is in balance with the amount and quality of the experimental data. If this is the case the model is said to be identifiable for the given data. If a model turns out to be unidentifiable, two steps can be taken. Either additional experiments need to be conducted, or the model has to be simplified.     

Sensitivity functions in the estimation of parameters of cellular exchange

The analysis of multiple-indicator dilution curves to estimate the rates of transport of ions and substrates across the sarcolemma of myocardial cells requires the formulation of models for the blood-interstitial fluid-cell exchanges. The fitting of models to the sets of experimental data is dependent on acquiring a large enough data set, in one physiological state, that there is at least as much information in the data as there are unknown model parameters to be determined. Inasmuch as data are necessarily noisy, redundancy of data and overdetermination of the unknowns are highly desirable. Sensitivity functions are useful in demonstrating which portions of the data relate to which unknown parameters. They are also useful in adjusting model parameters to fit the model to the data, and therefore in parameter evaluation.     

Adjoint multi-start-based estimation of cardiac hyperelastic material parameters using shear data

Cardiac muscle tissue during relaxation is commonly modeled as a hyperelastic material with strongly nonlinear and anisotropic stress response. Adapting the behavior of such a model to experimental or patient data gives rise to a parameter estimation problem which involves a significant number of parameters. Gradient-based optimization algorithms provide a way to solve such nonlinear parameter estimation problems with relatively few iterations, but require the gradient of the objective functional with respect to the model parameters. This gradient has traditionally been obtained using finite differences, the calculation of which scales linearly with the number of model parameters, and introduces a differencing error. By using an automatically derived adjoint equation, we are able to calculate this gradient more efficiently, and with minimal implementation effort. We test this adjoint framework on a least squares fitting problem involving data from simple shear tests on cardiac tissue samples. A second challenge which arises in gradient-based optimization is the dependency of the algorithm on a suitable initial guess. We show how a multi-start procedure can alleviate this dependency. Finally, we provide estimates for the material parameters of the Holzapfel and Ogden strain energy law using finite element models together with experimental shear data.     

Inverse analysis of constitutive models: biological soft tissues

The paper describes a procedure for estimating the material parameters of biological soft tissue by fitting model prediction to experimental load-deformation data. This procedure minimizes the error between data and theoretical model prediction through systematically adjusting the parameters in the latter. The procedure uses commercially available software and is not specific to any particular model; nevertheless, for illustration purposes, we employ a six parameter fibril-reinforced poroelastic cartilage model. We are able to estimate any and all of these parameters by the procedure. Convergence of the parameters and convergence of the arbitrary initial stress relaxation to the data was demonstrated in all cases. Though we illustrate the optimization procedure here for unconfined compression only, it can be adapted easily to other experimental configurations such as confined compression, indentation and tensile test. Furthermore, the procedure can be applied in other areas of biomechanics where material parameters need to be extracted from experimental data.     

Stochastic analogues of deterministic single-species population models

Although single-species deterministic difference equations have long been used in modeling the dynamics of animal populations, little attention has been paid to how stochasticity should be incorporated into these models. By deriving stochastic analogues to difference equations from first principles, we show that the form of these models depends on whether noise in the population process is demographic or environmental. When noise is demographic, we argue that variance around the expectation is proportional to the expectation. When noise is environmental the variance depends in a non-trivial way on how variation enters into model parameters, but we argue that if the environment affects the population multiplicatively then variance is proportional to the square of the expectation. We compare various stochastic analogues of the Ricker map model by fitting them, using maximum likelihood estimation, to data generated from an individual-based model and the weevil data of Utida. Our demographic models are significantly better than our environmental models at fitting noise generated by population processes where noise is mainly demographic. However, the traditionally chosen stochastic analogues to deterministic models--additive normally distributed noise and multiplicative lognormally distributed noise--generally fit all data sets well. Thus, the form of the variance does play a role in the fitting of models to ecological time series, but may not be important in practice as first supposed.     

Parameter redundancy in discrete state‐space and integrated models

Discrete state‐space models are used in ecology to describe the dynamics of wild animal populations, with parameters, such as the probability of survival, being of ecological interest. For a particular parametrization of a model it is not always clear which parameters can be estimated. This inability to estimate all parameters is known as parameter redundancy or a model is described as nonidentifiable. In this paper we develop methods that can be used to detect parameter redundancy in discrete state‐space models. An exhaustive summary is a combination of parameters that fully specify a model. To use general methods for detecting parameter redundancy a suitable exhaustive summary is required. This paper proposes two methods for the derivation of an exhaustive summary for discrete state‐space models using discrete analogues of methods for continuous state‐space models. We also demonstrate that combining multiple data sets, through the use of an integrated population model, may result in a model in which all parameters are estimable, even though models fitted to the separate data sets may be parameter redundant.     

Is Model Fitting Necessary for Model-Based fMRI?

Model-based analysis of fMRI data is an important tool for investigating the computational role of different brain regions. With this method, theoretical models of behavior can be leveraged to find the brain structures underlying variables from specific algorithms, such as prediction errors in reinforcement learning. One potential weakness with this approach is that models often have free parameters and thus the results of the analysis may depend on how these free parameters are set. In this work we asked whether this hypothetical weakness is a problem in practice. We first developed general closed-form expressions for the relationship between results of fMRI analyses using different regressors, e.g., one corresponding to the true process underlying the measured data and one a model-derived approximation of the true generative regressor. Then, as a specific test case, we examined the sensitivity of model-based fMRI to the learning rate parameter in reinforcement learning, both in theory and in two previously-published datasets. We found that even gross errors in the learning rate lead to only minute changes in the neural results. Our findings thus suggest that precise model fitting is not always necessary for model-based fMRI. They also highlight the difficulty in using fMRI data for arbitrating between different models or model parameters. While these specific results pertain only to the effect of learning rate in simple reinforcement learning models, we provide a template for testing for effects of different parameters in other models.     

Modelling avian growth with the Unified‐Richards: as exemplified by wader‐chick growth

Postnatal growth in birds is traditionally modelled by fitting three‐parameter models, namely the logistic, the Gompertz, or the von Bertalanffy models. The purpose of this paper is to address the utility of the Unified‐Richards (U‐Richards) model. We draw attention to two forms of the U‐Richards and lay down a set of recommendations for the analysis of bird growth, in order to make this model and the methods more accessible. We examine the behaviour of the four parameters in each model form and the four derived measurements, and we show that all are easy to interpret, and that each parameter controls a single curve characteristic. The two parameters that control the inflection point, enable us to compare its placement in two dimensions, 1) inflection value (mass or length at inflection) and 2) inflection time (time since hatching), between data sets (e.g. between biometrics or between species). We also show how the parameter controlling growth rate directly presents us with the relative growth rate at inflection, and we demonstrate how one can compare growth rates across data sets. The three traditional models, where the inflection value is fixed (to a specific percentage of the upper asymptote), provide incompatible growth‐rate coefficients. One of the two forms of the U‐Richards model makes it possible to fix not only the upper asymptote (adult value), but also the intersection with the y‐axis (hatching value). Fitting the new model forms to data validates the usefulness of interpreting the inflection placement in addition to the growth rate. It also illustrated the advantages and limitations of constraining the upper asymptote (adult value) and the y‐axis intersection (hatching value) to fixed values. We show that the U‐Richards model can successfully replace some of the commonly used growth models, and we advocate replacing these with the U‐Richards when modelling bird growth.     

Moment fitting for parameter inference in repeatedly and partially observed stochastic biological models

The inference of reaction rate parameters in biochemical network models from time series concentration data is a central task in computational systems biology. Under the assumption of well mixed conditions the network dynamics are typically described by the chemical master equation, the Fokker Planck equation, the linear noise approximation or the macroscopic rate equation. The inverse problem of estimating the parameters of the underlying network model can be approached in deterministic and stochastic ways, and available methods often compare individual or mean concentration traces obtained from experiments with theoretical model predictions when maximizing likelihoods, minimizing regularized least squares functionals, approximating posterior distributions or sequentially processing the data. In this article we assume that the biological reaction network can be observed at least partially and repeatedly over time such that sample moments of species molecule numbers for various time points can be calculated from the data. Based on the chemical master equation we furthermore derive closed systems of parameter dependent nonlinear ordinary differential equations that predict the time evolution of the statistical moments. For inferring the reaction rate parameters we suggest to not only compare the sample mean with the theoretical mean prediction but also to take the residual of higher order moments explicitly into account. Cost functions that involve residuals of higher order moments may form landscapes in the parameter space that have more pronounced curvatures at the minimizer and hence may weaken or even overcome parameter sloppiness and uncertainty. As a consequence both deterministic and stochastic parameter inference algorithms may be improved with respect to accuracy and efficiency. We demonstrate the potential of moment fitting for parameter inference by means of illustrative stochastic biological models from the literature and address topics for future research.     

Are Assumptions about the Model Type Necessary in Reaction-Diffusion Modeling? A FRAP Application

At present, fluorescence recovery after photobleaching (FRAP) data are interpreted using various types of reaction-diffusion (RD) models: the model type is usually fixed first, and corresponding model parameters are inferred subsequently. In this article, we describe what we believe to be a novel approach for RD modeling without using any assumptions of model type or parameters. To the best of our knowledge, this is the first attempt to address both model-type and parameter uncertainties in inverting FRAP data. We start from the most general RD model, which accounts for a flexible number of molecular fractions, all mobile, with different diffusion coefficients. The maximal number of possible binding partners is identified and optimal parameter sets for these models are determined in a global search of the parameter-space using the Simulated Annealing strategy. The numerical performance of the described techniques was assessed using artificial and experimental FRAP data. Our general RD model outperformed the standard RD models used previously in modeling FRAP measurements and showed that intracellular molecular mobility can only be described adequately by allowing for multiple RD processes. Therefore, it is important to search not only for the optimal parameter set but also for the optimal model type.     

Evaluating automated parameter constraining procedures of neuron models by experimental and surrogate data

Neuron models, in particular conductance-based compartmental models, often have numerous parameters that cannot be directly determined experimentally and must be constrained by an optimization procedure. A common practice in evaluating the utility of such procedures is using a previously developed model to generate surrogate data (e.g., traces of spikes following step current pulses) and then challenging the algorithm to recover the original parameters (e.g., the value of maximal ion channel conductances) that were used to generate the data. In this fashion, the success or failure of the model fitting procedure to find the original parameters can be easily determined. Here we show that some model fitting procedures that provide an excellent fit in the case of such model-to-model comparisons provide ill-balanced results when applied to experimental data. The main reason is that surrogate and experimental data test different aspects of the algorithm’s function. When considering model-generated surrogate data, the algorithm is required to locate a perfect solution that is known to exist. In contrast, when considering experimental target data, there is no guarantee that a perfect solution is part of the search space. In this case, the optimization procedure must rank all imperfect approximations and ultimately select the best approximation. This aspect is not tested at all when considering surrogate data since at least one perfect solution is known to exist (the original parameters) making all approximations unnecessary. Furthermore, we demonstrate that distance functions based on extracting a set of features from the target data (such as time-to-first-spike, spike width, spike frequency, etc.)—rather than using the original data (e.g., the whole spike trace) as the target for fitting—are capable of finding imperfect solutions that are good approximations of the experimental data.     

A Probabilistic Model for Indel Evolution: Differentiating Insertions from Deletions

Insertions and deletions (indels) are common molecular evolutionary events. However, probabilistic models for indel evolution are under-developed due to their computational complexity. Here, we introduce several improvements to indel modeling: 1) While previous models for indel evolution assumed that the rates and length distributions of insertions and deletions are equal, here we propose a richer model that explicitly distinguishes between the two; 2) we introduce numerous summary statistics that allow approximate Bayesian computation-based parameter estimation; 3) we develop a method to correct for biases introduced by alignment programs, when inferring indel parameters from empirical data sets; and 4) using a model-selection scheme, we test whether the richer model better fits biological data compared with the simpler model. Our analyses suggest that both our inference scheme and the model-selection procedure achieve high accuracy on simulated data. We further demonstrate that our proposed richer model better fits a large number of empirical data sets and that, for the majority of these data sets, the deletion rate is higher than the insertion rate.