Analytical and Simulation Results for Stochastic Fitzhugh-Nagumo Neurons and Neural Networks

An analytical approach is presented for determining the response of a neuron or of the activity in a network of connected neurons, represented by systems of nonlinear ordinary stochastic differential equations—the Fitzhugh-Nagumo system with Gaussian white noise current. For a single neuron, five equations hold for the first- and second-order central moments of the voltage and recovery variables. From this system we obtain, under certain assumptions, five differential equations for the means, variances, and covariance of the two components. One may use these quantities to estimate the probability that a neuron is emitting an action potential at any given time. The differential equations are solved by numerical methods. We also perform simulations on the stochastic Fitzugh-Nagumo system and compare the results with those obtained from the differential equations for both sustained and intermittent deterministic current inputs withsuperimposed noise. For intermittent currents, which mimic synaptic input, the agreement between the analytical and simulation results for the moments is excellent. For sustained input, the analytical approximations perform well for small noise as there is excellent agreement for the moments. In addition, the probability that a neuron is spiking as obtained from the empirical distribution of the potential in the simulations gives a result almost identical to that obtained using the analytical approach. However, when there is sustained large-amplitude noise, the analytical method is only accurate for short time intervals. Using the simulation method, we study the distribution of the interspike interval directly from simulated sample paths. We confirm that noise extends the range of input currents over which (nonperiodic) spike trains may exist and investigate the dependence of such firing on the magnitude of the mean input current and the noise amplitude. For networks we find the differential equations for the means, variances, and covariances of the voltage and recovery variables and show how solving them leads to an expression for the probability that a given neuron, or given set of neurons, is firing at time t. Using such expressions one may implement dynamical rules for changing synaptic strengths directly without sampling. The present analytical method applies equally well to temporally nonhomogeneous input currents and is expected to be useful for computational studies of information processing in various nervous system centers.     

Determination of active concentrations and association and dissociation rate constants of interacting biomolecules: an analytical solution to the theory for kinetic and mass transport limitations in biosensor technology and its experimental verification

Accurate determination of kinetic rate constants for interacting biomolecules requires knowledge of the active concentrations of the participating molecules. Also, in other biomedical and clinical applications, sensitive, precise and accurate methods are needed to determine the concentration of biologically active molecules, which frequently constitute only a fraction of the total molecular pool. Here we report a novel development of the approach to determining active concentrations based on surface plasmon resonance (SPR) technology. The method relies on changes in binding rates with varying flow rates under conditions of partial mass transport, and does not require standards of known concentrations, given that the molecular mass of the molecule of interest is known. We introduce an analytical solution to the differential equations describing the formation of a 1:1 bimolecular complex, taking into account both the association and dissociation reactions, under partial mass transport limitations. This solution can be used in global fitting to binding curves obtained at different flow rates. The accuracy, precision, and sensitivity of this approach were determined in experiments involving binding of tyrosine-phosphorylated recombinant proteins to anti-phosphotyrosine antibodies, where the active concentration could be determined independently by in vitro phosphorylation with (33)P. There was an excellent agreement between the active concentrations determined by the analytical SPR-based method and by determination of the level of radioactivity of the phosphorylated protein. The SPR-based method allows determination of protein concentrations at picomolar levels. A procedure for accurate determinations of association and dissociation rate constants, based on the analytical solution of the mass transport and binding theory, is outlined.     

On the Probability of Random Genetic Mutations for Various Types of Tumor Growth

In this work, we consider the problem of estimating the probability for a specific random genetic mutation to be present in a tumor of a given size. Previous mathematical models have been based on stochastic methods where the tumor was assumed to be homogeneous and, on average, growing exponentially. In contrast, we are able to obtain analytical results for cases where the exponential growth of cancer has been replaced by other, arguably more realistic types of growth of a heterogeneous tumor cell population. Our main result is that the probability that a given random mutation will be present by the time a tumor reaches a certain size, is independent of the type of curve assumed for the average growth of the tumor, at least for a general class of growth curves. The same is true for the related estimate of the expected number of mutants present in a tumor of a given size, if mutants are indeed present.     

Geometry of mediating protein affects the probability of loop formation in DNA

Recent single molecule experiments have determined the probability of loop formation in DNA as a function of the DNA contour length for different types of looping proteins. The optimal contour length for loop formation as well as the probability density functions have been found to be strongly dependent on the type of looping protein used. We show, using Monte Carlo simulations and analytical calculations, that these observations can be replicated using the wormlike-chain model for double-stranded DNA if we account for the nonzero size of the looping protein. The simulations have been performed in two dimensions so that bending is the only mode of deformation available to the DNA while the geometry of the looping protein enters through a single variable which is representative of its size. We observe two important effects that seem to directly depend on the size of the enzyme: 1), the overall propensity of loop formation at any given value of the DNA contour length increases with the size of the enzyme; and 2), the contour length corresponding to the first peak as well as the first well in the probability density functions increases with the size of the enzyme. Additionally, the eigenmodes of the fluctuating shape of the looped DNA calculated from simulations and theory are in excellent agreement, and reveal that most of the fluctuations in the DNA occur in regions of low curvature.     

Probability distributions of ancestries and genealogical distances on stochastically generated rooted binary trees

The stationary birth-only, or Yule-Furry, process for rooted binary trees has been analysed with a view to developing explicit expressions for two fundamental statistical distributions: the probability that a randomly selected leaf is preceded by N nodes, or “ancestors”, and the probability that two randomly selected leaves are separated by N nodes. For continuous-time Yule processes, the first of these distributions is presented in closed analytical form as a function of time, with time being measured with respect to the moment of “birth” of the common ancestor (which is essentially inaccessible to phylogenetic analysis), or with respect to the instant at which the first bifurcation occurred.The second distribution is shown to follow in an iterative manner from a hierarchy of second-order ordinary differential equations.For Yule trees of a given number n of tips, expressions have been derived for the mean and variance for each of these distributions as functions of n, as well as for the distributions themselves.In addition, it is shown how the methods developed to obtain these distributions can be employed to find, with minor effort, expressions for the expectation values of two statistics on Yule trees, the Sackin index (sum over all root-to-leaf distances), and the sum over all leaf-to-leaf distances.     

Dual polarisation interferometry characterisation of DNA immobilisation and hybridisation detection on a silanised support

Dual polarisation interferometry is an analytical technique that allows the simultaneous determination of thickness, density and mass of a biological layer on a sensing waveguide surface in real time. We evaluated, for the first time, the ability of this technique to characterise the covalent immobilisation of single stranded probe DNA and the selective detection of target DNA hybridisation on a silanised support. Two immobilisation strategies have been evaluated: direct attachment of the probe molecule and a more complex chemistry employing a 1,2 homobifunctional crosslinker molecule. With this technique we demonstrate it was possible to determine probe orientation and measure probe coverage at different stages of the immobilisation process in real time and in a single experiment. In addition, by measuring simultaneously changes in thickness and density of the probe layer upon hybridisation of target DNA, it was possible to directly elucidate the impact that probe mobility had on hybridisation efficiency. Direct covalent attachment of an amine modified 19 mer resulted in a thickness change of 0.68 nm that was consistent with multipoint attachment of the probe molecule to the surface. Blocking with BSA formed a dense layer of protein molecules that absorbed between the probe molecules on the surface. The observed hybridisation efficiency to target DNA was approximately 35%. No further significant reorientation of the probe molecule occurred upon hybridisation. The initial thickness of the probe layer upon attachment to the crosslinker molecule was 0.5 nm. Significant reorientation of the probe molecule surface normal occurred upon hybridisation to target DNA. This indicated that the probe molecule had greater mobility to hybridise to target DNA. The observed hybridisation efficiency for target DNA was approximately 85%. The results show that a probe molecule attached to the surface via a crosslinker group is better able to hybridise to target DNA due to its greater mobility.     

Membrane transport models with fast and slow reactions: General analytical solution for a single relaxation

Summary Membrane transport models are usually expressed on the basis of chemical kinetics. The states of a transporter are related by rate constants, and the time-dependent changes of these states are given by linear differential equations of first order. To calculate the time-dependent transport equation, it is necessary to solve a system of differential equations which does not have a general analytical solution if there are more than five states. Since transport measurements in a complex system rarely provide all the time constants because some of them are too rapid, it is more appropriate to obtain approximate analytical solutions, assuming that there are fast and slow reaction steps. The states of the fast steps are related by equilibrium constants, thus permitting the elimination of their differential equations and leaving only those for the slow steps. With a system having only two slow steps, a single differential equation is obtained and the state equations have a single relaxation. Initial conditions for the slow reactions are determined after the perturbation which redistribute the states related by fast reactions. Current and zero-trans uptake equations are calculated. Curve fitting programs can be used to implement the general procedure and obtain the model parameters.     

Stochastic model of population growth and spread

A stochastic model for the population regulated by logistic growth and spreading in a given region of two-or three-dimensional space has been introduced. For many-species population the interactions among the species have also been icorporated in this model. From the random variables that describe stochastic processes of a Wiener type the space-dependent random population densities have been formed and shown to satisfy the Langevin equations. The Fokker-Planck equation corresponding to these Langevin equations has been approximately solved for the transition probability of the population spreading and it has been found that such approximate expressions of the transition probability depend on the solutions of the deterministic equations of the diffusion model with logistic growth and interactions. Also, the stationary or equilibrium solutions of the Fokker-Planck equation together with the special discussion on the pattern of single-species population spreading have been made.     

New mathematical model to estimate tissue blood perfusion, thermal contact resistance and core temperature

Analytical solutions were developed based on the Green's function method to describe heat transfer in tissue including the effects of blood perfusion. These one-dimensional transient solutions were used with a simple parameter estimation technique and experimental measurements of temperature and heat flux at the surface of simulated tissue. It was demonstrated how such surface measurements can be used during step changes in the surface thermal conditions to estimate the value of three important parameters: blood perfusion (w(b)), thermal contact resistance (R"), and core temperature of the tissue (T(core)). The new models were tested against finite-difference solutions of thermal events on the surface to show the validity of the analytical solution. Simulated data was used to demonstrate the response of the model in predicting optimal parameters from noisy temperature and heat flux measurements. Finally, the analytical model and simple parameter estimation routine were used with actual experimental data from perfusion in phantom tissue. The model was shown to provide a very good match with the data curves. This demonstrated the first time that all three of these important parameters (w(b), R", and T(core)) have simultaneously been estimated from a single set of thermal measurements at the surface of tissue.     

Polymerase chain reaction: a Markov process approach

A probabilistic approach to the kinetics of the polymerase chain reaction (PCR) is developed. The approach treats the primer extension step of PCR as a microscopic Markov process in which the molecules of deoxy-nucleoside triphosphate (dNTP) are bound to the 3' end of the primer strand one at a time. The binding probability rates are prescribed by combinatorial rules in accord with the microscopic chemical kinetics. As an example, a simple model based on this approach is proposed and analysed, and an exact solution for the probability distribution of lengths of synthesized DNA strands is found by analytical means. Using this solution, it is demonstrated that the model is able to reproduce the main features of PCR, such as extreme sensitivity to the variation of control parameters and the existence of an amplification plateau. A multidimensional optimization technique is used to find numerically the optimum values of control parameters which maximize the yield of the target sequence for a given PCR run while minimizing the overall run time.     

Probability distribution of haplotype frequencies under the two-locus Wright-Fisher model by diffusion approximation

The probability distribution of haplotype frequencies in a population, and the way it is influenced by genetical forces such as recombination, selection, random drift ...is a question of fundamental interest in population genetics. For large populations, the distribution of haplotype frequencies for two linked loci under the classical Wright-Fisher model is almost impossible to compute because of numerical reasons. However the Wright-Fisher process can in such cases be approximated by a diffusion process and the transition density can then be deduced from the Kolmogorov equations. As no exact solution has been found for these equations, we developed a numerical method based on finite differences to solve them. It applies to transient states and models including selection or mutations. We show by several tests that this method is accurate for computing the conditional joint density of haplotype frequencies given that no haplotype has been lost. We also prove that it is far less time consuming than other methods such as Monte Carlo simulations.     

Coupling Between Folding and Ionization Equilibria: Effects of pH on the Conformational Preferences of Polypeptides

A new approach to the conformational study of polypeptides is presented. It considers explicitly the coupling between the conformation of the molecule and the ionization equilibria at a given pH value. Calculations of the solvation free energy and free energy of ionization of a 17-residue polypeptide are carried out using a fast multigrid boundary element method (MBE). The MBE method uses an adaptive tessellation of the molecular surface by boundary elements with non-regular size to solve the Poisson equation rapidly, and with a high degree of accuracy. The MBE method is integrated into the ECEPP (Empirical Conformational Energy Program for Peptides) algorithm to compute the coupling between the ionization state and the conformation of the molecule.This approach has been applied to study the conformational preference of a short polypeptide for which the available NMR and CD experimental data indicate that conformations containing a right-handed α-helical segment are energetically more favorable at low values of pH. The results of calculations using the present method agree quite well with experiments, in contrast to previous applications with standard techniques (using pre-assigned charges at each pH) that were not able to reproduce the experimental findings. Also, it is shown how the coupling to the conformation leads to different degrees of ionization of a given type of residue, for example glutamic acid, at different positions in the amino acid sequence, at any given pH. The results of this study provide a sound basis to discuss the origin of the stability of polypeptide conformations, and its dependence on the environmental conditions.     

Comparison of analytically and numerically derived hydrogen exchange-rate distribution functions

Hydrogen exchange-rate probability density functions for lysozyme have been derived by numerical Laplace inversion with the computer program CONTIN. The resulting solution set includes a smooth bimodal solution in agreement with previous analytical results together with a smooth three-peak solution. Numerical analysis of lysozyme hydrogen-exchange data in glycerol/water cosolvent mixtures confirms the previous assignment of the slow-exchange peak to an exchange mechanism involving reversible unfolding. Physicochemical constrations that can reduce the size of the solution set are described. The results are compared with those obtained from previous analytical methods and the limitations of the discrete class and analytical appraches are discussed.     

Hydrodynamics of bolus flow—An analytical approach to blood flow in capillaries

A model which was used by Prothero and Burton to simulate a particular configuration in capillary blood vessels is investigated from a hydrodynamic point of view. In this model, the erythrocytes are approximated by rigid pistons, and plasma is assumed to be an incompressible Newtonian fluid. An order of magnitude analysis using the physiologically realistic values for various parameters reduces the exact equations of motion to an equation describing the creeping motion of the fluid. An analytical approach to the solution of the equation is proposed and some results are reported here. The solution of the flow field is given in terms of a stream function which is represented by two infinite series composed of known functions. Two coupled infinite systems of algebraic equations determining the coefficients of the two series have been derived. This method of solution is proposed as an alternative to the entirely numerical procedure of solving the similar problem proposed by Bugliarelloet al. A limiting case of large aspect ratio (the ratio of the axial spacing of the two successive erythrocytes to the capillary diameter) is studied and the solution, valid away from the erythrocyte surface, has been obtained in simple form. It resembles the classical Poisenille flow, but the pressure gradient is related to the erythrocyte speed.     

NMR spin exchange kinetics at equilibrium in membrane transport and enzyme systems

A set of differential equations is formulated to describe the rapid exchange (time scale, approximately 0.01 to approximately 10 s) of a labelled solute across the membranes of cells in suspension. The labelling is achieved with nuclear magnetic resonance by exposure of the system to a high intensity radio-frequency pulse, and the excited nuclei relax to the equilibrium state with a short half life. An analytical expression for the decay of the magnetic resonance signal is presented; the solution involves the determination of eigenvalues, of an array of Laplace-Carson transformed differential equations, by use of the general solution of a quartic polynomial. Simulations of the behaviour of the exchange system using various conditions of cell number, rate constants and nuclear magnetic relaxation times are presented. The marked concentration dependence of the extent of reaction at a given time has not previously been reported for nuclear magnetic resonance exchange systems and is a feature anticipated from the known saturability of several membrane transport systems including glucose transport into human erythrocytes. The theory is readily generalized to other model systems by appropriate reinterpretation of the physical meaning of various parameters; the general form of the solution holds in many biological contexts other than membrane transport and includes equilibrium enzyme kinetics.     

Test of the Inverse Monte Carlo Method for the Calculation of Interatomic Potential Energies in Atomic Liquids

Abstract

The principle purpose of this paper is to demonstrate the use of the Inverse Monte Carlo technique for calculating pair interaction energies in monoatomic liquids from a given equilibrium property. This method is based on the mathematical relation between transition probability and pair potential given by the fundamental equation of the “importance sampling” Monte Carlo method. In order to have well defined conditions for the test of the Inverse Monte Carlo method a Metropolis Monte Carlo simulation of a Lennard Jones liquid is carried out to give the equilibrium pair correlation function determined by the assumed potential. Because an equilibrium configuration is prerequisite for an Inverse Monte Carlo simulation a model system is generated reproducing the pair correlation function, which has been calculated by the Metropolis Monte Carlo simulation and therefore representing the system in thermal equilibrium. This configuration is used to simulate virtual atom displacements. The resulting changes in atom distribution for each single simulation step are inserted in a set of non-linear equations defining the transition probability for the virtual change of configuration. The solution of the set of equations for pair interaction energies yields the Lennard Jones potential by which the equilibrium configuration has been determined.     

Adhesion frequency assay for in situ kinetics analysis of cross-junctional molecular interactions at the cell-cell interface

The micropipette adhesion assay was developed in 1998 to measure two-dimensional (2D) receptor-ligand binding kinetics. The assay uses a human red blood cell (RBC) as adhesion sensor and presenting cell for one of the interacting molecules. It employs micromanipulation to bring the RBC into contact with another cell that expresses the other interacting molecule with precisely controlled area and time to enable bond formation. The adhesion event is detected as RBC elongation upon pulling the two cells apart. By controlling the density of the ligands immobilized on the RBC surface, the probability of adhesion is kept in mid-range between 0 and 1. The adhesion probability is estimated from the frequency of adhesion events in a sequence of repeated contact cycles between the two cells for a given contact time. Varying the contact time generates a binding curve. Fitting a probabilistic model for receptor-ligand reaction kinetics to the binding curve returns the 2D affinity and off-rate. The assay has been validated using interactions of Fcγ receptors with IgG Fc, selectins with glycoconjugate ligands, integrins with ligands, homotypical cadherin binding, T cell receptor and coreceptor with peptide-major histocompatibility complexes. The method has been used to quantify regulations of 2D kinetics by biophysical factors, such as the membrane microtopology, membrane anchor, molecular orientation and length, carrier stiffness, curvature, and impingement force, as well as biochemical factors, such as modulators of the cytoskeleton and membrane microenvironment where the interacting molecules reside and the surface organization of these molecules. The method has also been used to study the concurrent binding of dual receptor-ligand species, and trimolecular interactions using a modified model. The major advantage of the method is that it allows study of receptors in their native membrane environment. The results could be very different from those obtained using purified receptors. It also allows study of the receptor-ligand interactions in a sub-second timescale with temporal resolution well beyond the typical biochemical methods. To illustrate the micropipette adhesion frequency method, we show kinetics measurement of intercellular adhesion molecule 1 (ICAM-1) functionalized on RBCs binding to integrin α(L)β(2) on neutrophils with dimeric E-selectin in the solution to activate α(L)β(2).