A fast and reliable computational method for estimating population genetic parameters

The estimation of ancestral and current effective population sizes in expanding populations is a fundamental problem in population genetics. Recently it has become possible to scan entire genomes of several individuals within a population. These genomic data sets can be used to estimate basic population parameters such as the effective population size and population growth rate. Full-data-likelihood methods potentially offer a powerful statistical framework for inferring population genetic parameters. However, for large data sets, computationally intensive methods based upon full-likelihood estimates may encounter difficulties. First, the computational method may be prohibitively slow or difficult to implement for large data. Second, estimation bias may markedly affect the accuracy and reliability of parameter estimates, as suggested from past work on coalescent methods. To address these problems, a fast and computationally efficient least-squares method for estimating population parameters from genomic data is presented here. Instead of modeling genomic data using a full likelihood, this new approach uses an analogous function, in which the full data are replaced with a vector of summary statistics. Furthermore, these least-squares estimators may show significantly less estimation bias for growth rate and genetic diversity than a corresponding maximum-likelihood estimator for the same coalescent process. The least-squares statistics also scale up to genome-sized data sets with many nucleotides and loci. These results demonstrate that least-squares statistics will likely prove useful for nonlinear parameter estimation when the underlying population genomic processes have complex evolutionary dynamics involving interactions between mutation, selection, demography, and recombination.     

Least-squares estimation of batch culture kinetic parameters

This article concerns the development of a simple and effective least-squares procedure for estimating the kinetic parameters in Monod expressions from batch culture data. The basic approach employed in this work was to translate the problem of parameter estimation to a mathematical model containing a single decision variable. The resulting model was then solved by an efficient one-dimensional search algorithm which can be adapted to any microcomputer or advanced programmable calculator. The procedure was tested on synthetic data (substrate concentrations) with different types and levels of error. The effect of endogeneous respiration on the estimated values of the kinetic parameters was also assessed. From the results of these analyses the least-squares procedure developed was concluded to be very effective.     

Parameter estimation and sensitivity analysis of a nonlinearly elastic static lung model

A model for the static pressure-volume behavior of the lung parenchyma based on a pseudo-elastic strain energy function was tested. Values of the model parameters and their variances were estimated by an optimal least-squares fit of the model-predicted pressures to the corresponding data from excised, saline-filled dog lungs. Although the model fit data from twelve lungs very well, the coefficients of variation for parameter values differed greatly. To analyze the sensitivity of the model output to its parameters, we examined an approximate Hessian, H, of the least-squares objective function. Based on the determinant and condition number of H, we were able to set formal criteria for choosing the most reliable estimates of parameter values and their variances. This in turn allowed us to specify a normal range of parameter values for these dog lungs. Thus the model not only describes static pressure-volume data, but also uses the data to estimate parameters from a fundamental constitutive equation. The optimal parameter estimation and sensitivity analysis developed here can be widely applied to other physiologic systems.     

Evaluation and propagation of confidence intervals in nonlinear, asymmetrical variance spaces. Analysis of ligand-binding data

Two problems that are often overlooked in studies employing nonlinear least-squares techniques for parameter estimation are confidence-interval estimation and propagation. When the parameters are correlated, the variance space and consequently the confidence intervals are nonlinear and asymmetrical. The presented mathematical method for the evaluation of confidence intervals and error propagation addresses these problems. The examples employed to demonstrate these methods include linear least-squares and the nonlinear least-squares analysis of ligand-binding problems, such as hormone receptor interactions and oxygen binding to human hemoglobin. The mathematical procedures have proven very useful for analyzing the molecular mechanism of cooperativity in human hemoglobin (Johnson, M. L., and G. K. Ackers, 1982. Biochemistry 21:201-211).     

Identification in an anaerobic batch system: global sensitivity analysis, multi-start strategy and optimization criterion selection

Several mathematical models have been developed in anaerobic digestion systems and a variety of methods have been used for parameter estimation and model validation. However, structural and parametric identifiability questions are relatively seldom addressed in the reported AD modeling studies. This paper presents a 3-step procedure for the reliable estimation of a set of kinetic and stoichiometric parameters in a simplified model of the anaerobic digestion process. This procedure includes the application of global sensitivity analysis, which allows to evaluate the interaction among the identified parameters, multi-start strategy that gives a picture of the possible local minima and the selection of optimization criteria or cost functions. This procedure is applied to the experimental data collected from a lab-scale sequencing batch reactor. Two kinetic parameters and two stoichiometric coefficients are estimated and their accuracy was also determined. The classical least-squares cost function appears to be the best choice in this case study.     

Power-law modeling based on least-squares minimization criteria.

The power-law formalism has been successfully used as a modeling tool in many applications. The resulting models, either as Generalized Mass Action or as S-systems models, allow one to characterize the target system and to simulate its dynamical behavior in response to external perturbations and parameter changes. The power-law formalism was first derived as a Taylor series approximation in logarithmic space for kinetic rate-laws. The especial characteristics of this approximation produce an extremely useful systemic representation that allows a complete system characterization. Furthermore, their parameters have a precise interpretation as local sensitivities of each of the individual processes and as rate-constants. This facilitates a qualitative discussion and a quantitative estimation of their possible values in relation to the kinetic properties. Following this interpretation, parameter estimation is also possible by relating the systemic behavior to the underlying processes. Without leaving the general formalism, in this paper we suggest deriving the power-law representation in an alternative way that uses least-squares minimization. The resulting power-law mimics the target rate-law in a wider range of concentration values than the classical power-law. Although the implications of this alternative approach remain to be established, our results show that the predicted steady-state using the least-squares power-law is closest to the actual steady-state of the target system.     

Methods for comparison of parameters from longitudinal rhythmometric models with multiple components

Multiple components linear least-squares methods have been proposed for the detection of periodic components in nonsinusoidal longitudinal time series. However, a proper test for comparison of parameters obtained from this method for two or more time series is not yet available. Accordingly, we propose two methods, one parametric and one nonparametric, to compare parameters from rhythmometric models with multiple components. The parametric method is based on techniques commonly and generally employed in linear regression analysis. The comparison of parameters among two or more time series is accomplished by the use of so-called dummy variables. The nonparametric method is based on bootstrap techniques. This approach basically tests if the difference in any given parameter obtained by fitting a model with the same periods to two different longitudinal time series differs from zero. This method calculates a confidence interval for the difference in the tested parameter. If this interval does not contain zero, it can be concluded that the parameters obtained from the two time series are different with high probability. An estimation of the p-value for the corresponding test can also be calculated. By the use of similar bootstrap techniques, confidence intervals can also be obtained for any parameter derived from the multiple component fit of several periods to nonsinusoidal longitudinal time series, including the orthophase (peak time), bathyphase (trough time), and global amplitude (difference between the maximum and the minimum) of the fitted model waveform. These methods represent a valuable tool for the comparison of rhythm parameters obtained by multiple component analysis, and they render this approach as a generally applicable one for waveform representation and detection of periodicities in nonsinusoidal, sparse, and noisy longitudinal time series sampled with either equidistant or unequidistant observations.     

HIV model parameter estimates from interruption trial data including drug efficacy and reservoir dynamics

Mathematical models based on ordinary differential equations (ODE) have had significant impact on understanding HIV disease dynamics and optimizing patient treatment. A model that characterizes the essential disease dynamics can be used for prediction only if the model parameters are identifiable from clinical data. Most previous parameter identification studies for HIV have used sparsely sampled data from the decay phase following the introduction of therapy. In this paper, model parameters are identified from frequently sampled viral-load data taken from ten patients enrolled in the previously published AutoVac HAART interruption study, providing between 69 and 114 viral load measurements from 3-5 phases of viral decay and rebound for each patient. This dataset is considerably larger than those used in previously published parameter estimation studies. Furthermore, the measurements come from two separate experimental conditions, which allows for the direct estimation of drug efficacy and reservoir contribution rates, two parameters that cannot be identified from decay-phase data alone. A Markov-Chain Monte-Carlo method is used to estimate the model parameter values, with initial estimates obtained using nonlinear least-squares methods. The posterior distributions of the parameter estimates are reported and compared for all patients.     

A user's experience with a standard non-linear regression program (BMDP3R).

A study was made to test and compare the behavior of a standard non-linear regression program (BMDP3R) in fitting data from six classical least-squares problems. The use of three program control parameters is discussed and four measures of regression failure are utilized to give a quantitative reference of success. Recommendations are given to aid the user of packaged programs in the parameter estimation of non-linear regression models.     

Time-dependent power spectral density estimation of surface electromyography during isometric muscle contraction: Methods and comparisons

This paper studies the time-dependent power spectral density (PSD) estimation of nonstationary surface electromyography (SEMG) signals and its application to fatigue analysis during isometric muscle contraction. The conventional time-dependent PSD estimation methods exhibit large variabilities in estimating the instantaneous SEMG parameters so that they often fail to identify the changing patterns of short-period SEMG signals and gauge the extent of fatigue in specific muscle groups. To address this problem, a time-varying autoregressive (TVAR) model is proposed in this paper to describe the SEMG signal, and then the recursive least-squares (RLS) and basis function expansion (BFE) methods are used to estimate the model coefficients and the time-dependent PSD. The instantaneous parameters extracted from the PSD estimation are evaluated and compared in terms of reliability, accuracy, and complexity. Experimental results on synthesized and real SEMG data show that the proposed TVAR-model-based PSD estimators can achieve more stable and precise instantaneous parameter estimation than conventional methods.     

Online estimation algorithm for a biaxial ankle kinematic model with configuration dependent joint axes

The kinematics of the human ankle is commonly modeled as a biaxial hinge joint model. However, significant variations in axis orientations have been found between different individuals and also between different foot configurations. For ankle rehabilitation robots, information regarding the ankle kinematic parameters can be used to estimate the ankle and subtalar joint displacements. This can in turn be used as auxiliary variables in adaptive control schemes to allow modification of the robot stiffness and damping parameters to reduce the forces applied at stiffer foot configurations. Due to the large variations observed in the ankle kinematic parameters, an online identification algorithm is required to provide estimates of the model parameters. An online parameter estimation routine based on the recursive least-squares (RLS) algorithm was therefore developed in this research. An extension of the conventional biaxial ankle kinematic model, which allows variation in axis orientations with different foot configurations had also been developed and utilized in the estimation algorithm. Simulation results showed that use of the extended model in the online algorithm is effective in capturing the foot orientation of a biaxial ankle model with variable joint axis orientations. Experimental results had also shown that a modified RLS algorithm that penalizes a deviation of model parameters from their nominal values can be used to obtain more realistic parameter estimates while maintaining a level of estimation accuracy comparable to that of the conventional RLS routine.     

Quantifying nonlinear anisotropic elastic material properties of biological tissue by use of membrane inflation

Determination of material parameters for soft tissue frequently involves regression of material parameters for nonlinear, anisotropic constitutive models against experimental data from heterogeneous tests. Here, parameter estimation based on membrane inflation is considered. A four parameter nonlinear, anisotropic hyperelastic strain energy function was used to model the material, in which the parameters are cast in terms of key response features. The experiment was simulated using finite element (FE) analysis in order to predict the experimental measurements of pressure versus profile strain. Material parameter regression was automated using inverse FE analysis; parameter values were updated by use of both local and global techniques, and the ability of these techniques to efficiently converge to a best case was examined. This approach provides a framework in which additional experimental data, including surface strain measurements or local structural information, may be incorporated in order to quantify heterogeneous nonlinear material properties.     

Determination of stoichiometric association constants by a non-iterative computational method

Computational methods for the estimation of stoichiometric associationconstants for multiple-ligand binding systems are currentlybased on non-linear least-squares regression analysis. Thesecomputational methods require sophisticated, iterative algorithmsto assure convergence to a solution, as well as initial parameterand error estimates. A simple procedure, called lambda-invariancetesting (LIT), was developed that provides a single-pass (non-iterative)estimation of stoichiometric association constants. The LITmethod was applied to simulated binding data containing Gaussianerror and to real data drawn from the literature. This methodprovided parameter estimates essentially equivalent to thoseobtained by least-squares regression analysis, with no initialparameter or error estimates required. Received on May 26, 1987; accepted on July 27, 1987     

A calculator program for clinical application of the Bayesian method of predicting plasma drug levels

A pharmacokinetic program that allows individualization of drug dosage regimens through the Bayesian method is described. The program, which is designed for the Hewlett-Packard HP-41 CV calculator, is based upon the one-compartment open model with either instantaneous or zero-order absorption. Individualized estimation of the patient's kinetic parameters (clearance and volume of distribution) is performed by analyzing the plasma levels measured in the patient as well as considering the population data of the drug. After estimating the individual kinetic parameters by the Bayesian method, the program predicts the dosage regimen that will elicit the desired peak and trough plasma levels at steady state. For comparison purposes, the least-squares estimates for clearance and volume of distribution are calculated, and dosage prediction can also be made on the basis of the least-squares estimates. The least-squares estimates can be used to calculate population pharmacokinetic parameters according to the Standard Two-Stage method. Several examples of clinical use of the program are presented. The examples refer to patients with classic hemophilia who were treated with Factor VIII concentrates. In these patients, the Bayesian kinetic parameters of Factor VIII have been estimated through the calculator program. The Bayesian parameter estimates generated by the HP-41 have been compared with those determined by a Bayesian program (ADVISE) designed for microcomputers.     

Quantification of two-dimensional NOE spectra via a combined linear and nonlinear least-squares fit

Summary Determining the volumes of peaks in 2D NMR spectra can be prohibitively difficult in cases of overlapping, broad lines. Deconvolution and parameter estimation can be attempted on either the time-domain or the frequency-domain data. We present a method of estimating spectral parameters from frequency-domain data, using a combination of Lorentzian and Gaussian lineshapes for reference lines. This approach combines a previously published method of projecting the data on a linear space spanned by reference lines with a nonlinear least-squares fitting algorithm. Comparison of this method with other published methods of frequency-domain deconvolution shows that it is both more precise and more accurate when estimating 2D volumes.     

Constrained least-squares optimization for robust estimation of center of rotation

This paper presents a new direct method for estimating the average center of rotation (CoR). An existing least-squares (LS) solution has been shown by previous works to have reduced accuracy for data with small range of motion (RoM). Alternative methods proposed to improve the CoR estimation use iterative algorithms. However, in this paper we show that with a carefully chosen normalization scheme, constrained least-squares solutions can perform as well as iterative approaches, even for challenging problems with significant noise and small RoM. In particular, enforcing the normalization constraint avoids poor fits near plane singularities that can affect the existing LS method. Our formulation has an exact solution, accounts for multiple markers simultaneously, and does not depend on manually-adjusted parameters. Simulation tests compare the method to four published CoR estimation techniques. The results show that the new approach has the accuracy of the iterative methods as well as the short computation time and repeatability of a least-squares solution. In addition, application of the new method to experimental motion capture data of the thumb carpometacarpal (CMC) joint yielded a more plausible CoR location compared to the previously reported LS solution and required less time than all four alternative techniques.     

A Method of Estimation of Change Points in Multiphasic Growth Models

The least-squares estimation (LSE) of change points and the phase parameters is considered for two growth models with piecewise continuously connected exponential and stationary phases. The calculation of a LSE is reduced to the solution of a finite number of simpler least-squares problems for which methods of calculation are well known.     

Parameter estimation for ligand binding systems kinetics applied to 1,25-dihydroxycholecalciferol

A mathematical analysis of the kinetics of the hormone-receptor interaction was applied to the 1,25-dihydroxycholecalciferol-intestinal receptor system. The exact analytical solution and the numerical integration of the kinetic equation were installed in a Statistical Analysis System (SAS) computer program to estimate the rate constants of the reaction. Estimates of the parameters obtained by these two methods are similar, demonstrating that the numerical integration can be combined with the nonlinear regression procedure for least-squares parameter fitting using a simple SAS program. This enables estimation of kinetics rate constants when the kinetic equation cannot be solved analytically. The ratio of the rate constants (ka/kd) found by the nonlinear procedure is close to the independently determined equilibrium (Scatchard) constant in the nonlinear analysis.     

Ionic-strength-dependent substrate inhibition of the lysis of Micrococcus luteus by hen egg-white lysozyme

The kinetics of lysis of Micrococcus luteus by hen egg-white lysozyme in dilute buffer media is characterized by pronounced substrate inhibition. This effect occurs within the complete pH range where lysozyme activity is detectable. The electrostatic potential of the negatively charged cell-wall proteoglycan increases with decreasing ionic strength, resulting in an enhanced affinity between proteoglycan and lysozyme and probably favouring multipoint substrate attachment. For the lysozyme-catalyzed hydrolysis of cell-wall proteoglycan three plausible mechanisms of substrate inhibition can be postulated. Two out of the three models fit our experimental data, the simplest of the two providing the most rigorous information on the kinetic parameters Km, V and Ki. Three graphical methods consistent with the chosen model were applied for preliminary parameter estimation and the constants obtained were compared to those from nonlinear least-squares analysis. If substrate inhibition is neglected it is shown that serious bias is imposed upon the parameters.     

Adaptive control of dissolved oxygen concentration in a bioreactor

A new adaptive DO (dissolved oxygen) concentration control algorithm considering DO electrode dynamics with response time delay has been developed. A system model with two time-varying parameters was used to relate the DO concentration with two control variables: air flow rate and agitation speed. Parameters of this model were estimated on-line using a regularized constant trace recursive least-squares method. An extended Kalman filter was used to remove the effect of noises from the DO concentration measurements and thus to improve control performance. A discrete one-step ahead control scheme was adopted to determine control actions based on the parameter estimation results. Experimental results showed that the new adaptive DO concentration control algorithm performed better than other algorithms tested, a PID controller and adaptive algorithms without the DO electrode dynamics.