The "inverse" problem is related to the determination of unknown causes on the bases of the observation of their effects. This is the opposite of the corresponding "direct" problem, which relates to the prediction of the effects generated by a complete description of some agencies. Typical inverse problems in electrocardiology are related to the modelling of the human heart functional structure from surface electrocardiogram signals (ECG) [1]; similar situations are encountered in magnetoencephalography (MEG) and electroencephalography (EEG) [2, 3]. In biology, a classical example of the "inverse" approach is the reconstruction of the three-dimensional structure of macromolecules, using either x-ray diffraction, nuclear magnetic resonance (NMR) or prediction models [4–6]. Another typical biological application of inverse approaches is the reconstruction of gene-regulatory networks [7, 8].

The solution of an inverse problem entails the construction of a mathematical model and takes the moves from a number of experimental data. In this respect, inverse problems are often ill-conditioned as the amount of experimental conditions available are often insufficient to unambiguously solve the mathematical model. Moreover, as model construction usually depends upon the minimization of specific functions, such as the system energy or the difference between the model prediction and some given experimental results, its solution does not necessarily lead to a single global optimal solution but to a set of optimal solutions, defining what is called the "Pareto optimal frontier" in the space of solutions [9]. Additional experimental constraints or theoretical methods are thus necessary to further select within the solutions set. Typical inverse problems concerns essentially the detailed determination of biochemical mechanisms underlying observed phenotypes, for example molecular abundances or morphological modifications.

In this work, we will attempt to solve an inverse problem which arose in the study of a signalling pathway. Compared to pathways of metabolic reactions, which are of a limited size comprising up to a few hundreds of proteins, signalling processes involve about 20% of the genome, i.e. thousands of expressed proteins [10], most still unidentified and of unknown function. Protein signalling networks spread information throughout the cell and mediate a number of fundamental processes [11–14]. The growing availability of reliable genomic and proteomic data, made it possible to build up protein interaction maps (PIMs) of increasing complexity. New high-throughput experimental and in silico technologies allow us to monitor protein-protein and genetic interactions: DNA and protein microarrays [15–17], two-hybrid systems [18–20], protein tagging techniques coupled with Mass Spectrometry [21, 22], phage display [23, 24]. In silico methods also allow us to describe protein-protein (p-p hereafter) interactions or the function of yet unclassified proteins: new p-p interactions might be found on the base of genomic sequence [25, 26], using data mining methodologies [27, 28], or predicting the composition of protein complexes [29]. In this respect it is worth mentioning a simple though successful method to detect new protein-protein interactions by a comparative genomic analysis of phylogenetic profiles: this approach is based on the assumption that interacting genes tend to co-evolve in different organisms [30, 31]. Protein's function can be predicted not only by sequence homology, but also on the basis of their relationships with other proteins whose role is already experimentally assessed [32, 33] or by orthology [34]. In order to model the time evolution of a signalling pathway it is necessary to know:

The mathematical form of the chemical interactions, the model parameters and even the network topology are often only partially known. This implies that model approximations and numerical estimates and, whenever possible, additional specific experimental measurements, are necessary to make a numerical simulation feasible and reliable. This is true whatever modelling techniques is used, such as differential equations [35, 36], cellular automata [37], Petri Nets [38] or other hybrid methods [39]. When creating a new model, before starting with numerical procedures, it is necessary to make a survey on all published kinetic data. These data may be found directly in the journal articles, which requires a thorough mining of the literature, or on in annotated databases, collecting and structuring information on p-p interactions.

Only at the end of this phase, further experimental activity and the techniques for parameter's estimate come into play: wherever possible, purposely designed experiments should be carried out in order to directly measure unknown kinetic parameters or to use these measures as constraints for the estimate's algorithm or to decide between alternative models. If new experiments cannot be done, the parameter estimate must rely just on literature data.

Protein interactions maps, partially stored in public databases, contain mainly qualitative information on the connectivity of intracellular p-p interactions, while quantitative data on the kinetics of interactions and reactions are still largely unavailable, except for enzyme kinetics. There are to date a number of public databases sites containing qualitative data on protein interaction maps:

It should be remarked that a great care has to be payed when dealing with qualitative data: they are often dependent on specific experimental conditions and most of them obtained in unicellular organisms. A straightforward extrapolation of these data to higher organisms is often quite unreliable [52]. Moreover, p-p interactions data in molecular networks are usually obtained from large scale or high-throughput experiments, where spurious interactions are very likely to be collected; computational validation techniques are thus needed to prune primary datasets [53, 54]. The same holds when one tries to translate genetic interactions into the corresponding p-p interactions: the two networks have quite different topological properties [55].

The situation is even worse when one analyzes quantitative p-p interactions data in public repositories: the total amount of experimentally-derived kinetic data is only a small percentage of what would be needed to characterize the topology data (i.e. the p-p interactions map). Furthermore, available kinetic constants are often extracted from a single publication where they were measured in vitro, while the kinetics of interactions is highly dependent on experimental set-up and environmental conditions, such as PH, temperature, concentration of other proteins in the cellular environment. It is always advisable to assume that the measured quantities indicate more realistically ranges rather than precise values and care must be used to insert these values into large-scale network models [56]. Nevertheless some investigation of biochemical reactions can anyway be carried out by taking into account the uncertainty of kinetic data [57, 58] and by using approximations where some values are missing [39].

This point, however, is already a major concern of the Systems Biology: several programs are being performed aimed at producing sets of validated data, homogeneously refered at specific organisms in well defined and standardized thermo-chemical conditions. The standardization of experimental data sets and of experimental models is the object of an intense debate in the Systems Biology community. There is a wide consensus on the need of standards but also on some drawbacks for a general use of standards as the best research framework in any case. Anyway the way towards a deeper and deeper though slow integration of existing datasets, modelling languages and methodologies appears to be set, as witnessed for example by the wider and wider use of SBML as a language to describe biochemical models, or by the integration of previously separated datasets into a single larger database compliant with new criteria established by international consortia. One example of the latter case is the HUPO – PSI initiative [59], aimed at establishing a common format to represent protein-protein interactions and to synchronize all the already existing databases, as it happened for the genome data: MINT, DIP, BIND and IntAct (see below) already implemented the PSI standard to publish molecular interactions.

p-p interactions in signalling pathways can be divided into two main categories: (a) binding interactions that involve no chemical modifications and (2) biochemical processing, such as phosphorylation and phosphatization. On one hand, the few public sources of kinetic data on binding protein interaction often provide only dissociation constants, i.e. values describing an equilibrium state that offer only partial information about the dynamics of the reaction. To our knowledge, only the KDBI database [60, 61] was specifically created to store binding and dissociation rate constants. Other repositories, such as MINT [62, 63] and BIND [64, 65] offer few examples of dissociation constants. On the other hand, biochemical modifications occur in enzymatic reactions, therefore kinetic data can be found in databases entirely devoted to enzymes, first of all Brenda [66, 67] where kinetic constants are specified for several organic substrates, and partially the above cited KDBI.

A further source of signalling pathways and of p-p interactions data, including the kinetic part, are the repositories of biochemical models, though in these models not all the kinetic parameters were measured experimentally and some of them had to be numerically estimated. Among them:

The measure of protein activation level is of paramount importance to monitor signalling processes. Several methods exist to quantitate the concentration of protein species, such as immunoblotting, ELISA, radioimmunoassay, protein arrays. If a cellular system is sampled several times over the duration of a given signalling process, a time series can be composed describing the time course of a concentration, for example that of a phosphorilated protein. Radioimmunoassays are very sensitive methods but are even complex, expensive and dangerous to set up; protein arrays offer the advantage of a high throughput approach, while ELISA and immunoblotting are easier to implement and, thus, widely used, though they allow a lower threshold of detection when a very low concentrations of radioactive compounds is present [79]. The experimental error of quantitative immunoblotting can be significantly reduced by computational techniques of data analysis, error estimate and simulation: these allow to monitor activated signalling pathways in real time and to discriminate between different models.

Enzymatic reactions can be monitored, nowadays, in a high throughput scale both in vivo and in vitro: this allows us to measure kinetic parameters characterizing fundamental steps in signalling pathways, such as binding and removal of phosphate groups by kinases and phosphatases. Bioreactors are widely used to perform enzymatic reactions and other biochemical processes but their use for a real time monitoring of products is limited by the sampling process. More recent modified reactors allow a real-time sampling of multiple reactions in vivo over a short reaction time: the reaction broth flows at constant velocity along a thin pipe where spilling at uniform space intervals corresponds to uniform time sampling. In this system the samples can be rapidly quenched and analyzed by mass-spectrometry techniques [80]. Also arrays of nanolitre wells can be used to follow the time course of multiple enzymatic process by the use of optical techniques such as fluorescence and bioluminescence [81]. The analysis of reaction mixtures by mass spectrometry methods makes the use of chromophores and radiolabelling unnecessary, since even the addition of a phosphate group to a large protein can be detected as a precise mass shift in the spectra. In vitro multiplexed assays can be performed on protein chips that are then directly analyzed by surface-enhanced laser desorption/ionization time-of-flight mass spectrometry (SELDI-TOF MS) to monitor enzyme activities [82]. Alternatively complex protein mixtures can be immobilized on micro-beads, where the enzymatic reactions can take place and be monitored by MALDI mass spectrometry [83]. A more difficult issue is to measure kinetic parameters describing binding of proteins without chemical processing, such as ligand-receptor interactions or the formation of protein complexes. Two techniques allow us to calculate kinetic rate constants of binding an unbinding by fitting measured response curves. The Surface Plasmon Resonance (SPR) allows us to measure kinetic constants in vitro in a label-free environment. One of the reactants is immobilized on the sensor surface usually coated with a thin gold film, while the other is free in solution: the behaviour of a polarized light beam hitting the surface in conditions of total internal reflection depends on the refractive index of the surface, that in turn depends on the binding state of the reactants. In essence the SPR measures the angle or the wavelength of the reflected light at which a resonance takes place between the light and the metal electrons: whose changes correspond to the amount of bound molecules. The SPR is already used for high-throughput measurements directly on protein arrays [84–86]. Using a completely different approach called fluorescence cross-correlation spectroscopy (FCS) the kinetics of binding can be quantified directly in living cells. Fluctuations of fluorescence signals can be detected in a very reduced volume, less than a femtolitre, and using a very low fluorophore concentration, up to 5 nM i.e. around 3 molecules/femtolitre, by the use of a tightly focused laser beam. The investigation of the autocorrelation function of the fluorescence signal provides information about the reaction kinetics, the diffusion rates and the equilibrium state. With FCS it is feasible to study at a single molecule level a ligand-receptor interaction with no need of any isotope labeling [87, 88].

When only a few kinetic parameters are available to implement a model of a signaling network, one might resort to attempting a "theoretical" estimate of these values. The attempt could be performed, in principle, by using an "inverse problem" approach, i.e. by optimizing the unknown parameters of a reaction's model in order to obtain the best possible agreement between simulated and experimental data.

This is the aim of the present work. We devise a methodological workflow (and the corresponding numerical and computational tools) to estimate the unknown reaction constants of a model signalling pathway by starting from (a) a given set of known data of reaction constants and (b) experimental results of the time course of some biochemical species involved in the reaction.

An intracellular signal transduction pathway in the neuronal cell was used as a model system to implement the proposed parameter's estimate procedure.

The chosen pathway is a protein network downstream of the neurotrophic receptors Trk and P75 [89], the Fas receptor regulating an apoptotic cascade [35], the EGF receptor expressed in the CNS [90, 91] and in PC12 cells [92, 93]. The network structure is based on current literature. The pathway can be divided into two main interconnected sub-systems: an apoptosis pathway and a neurotrophic receptors activated pathway. Neuronal apoptosis can be initiated in three different manners, all leading to the activation of executioner caspases, the effectors of the apoptotic process that kill the cell by irreversible proteolysis of critical cellular constituents: survival factor withdrawal, stress factors and receptor mediated signaling cascade [94, 35, 95]. In this model the survival factor withdrawal is taken into account by the connections between the two sub-networks, the apoptotic and the neurotrophic driven one, which includes the TRK, EGFR and P75NTR receptors the stress factors are also considered by the presence of a mitochondrion acting as a synthesis machinery for pro-apoptotic proteins (Fig. 1). The signaling pathways forming the network can be activated in several ways; in our model, we chose to trigger the signalling process by the activation of the receptors upstream of the pathways as a consequence of the binding of specific ligands.

The p-p interactions, such as molecular binding, phosphorylation/dephospshorylation or chemical transformations, are described using first order non-linear ordinary differential equations, which take into account also synthesis and degradation processes. The space variable is neglected in this model, since proteins are considered to be close enough to justify the approximation of a geometrical point. The release from the mitochondria was considered to be mathematically equivalent to an additional protein synthesis [94, 35]. In this model gene transcription was neglected, owing to the time scale chosen to simulate the temporal evolution of the system, within 60 minutes time. Reactions are treated as a one-step process. For binary activation and inactivation reactions, the following second order kinetics scheme was used, where protein A activates protein B:

The activation rate of protein B is : v _act = K _act [A][B]. In the case of binding reactions, resulting in the association/dissociation of protein complexes, the following one-step reaction scheme was used, resulting in a p ^th -order kinetics, where p equals the number of components of the protein complex C _i , with forward and reverse rate constants K and K^-1 respectively:

Thus the association rate is v _ass = K [C₁] [C₂]...[C _n ] and the dissociation rate is v _diss = K_-1[P].

Each of the N = 98 nodes of the network is described by the two independent variables P _i and x _i (i = l...N): the first refers to the total concentration of the protein species, the second to the concentration of the active fraction of that species. Each protein species i will thus follow a time evolution given by two coupled reactions:

where v_{prod, a}(v_cons,a), with a = x _i , P _i , represent production (consumption) reactions having the a-species as object. The complete system of equations describing the system assumes the following explicit mathematical structure:

N in the number of nodes, NP_i,jis the number of different interactions involving the nodes i and j, NC_i,j,ris the number of components when i is a protein in complex with protein j and the represent the different rate constants. The r index accounts for different interactions between nodes i and j, when existing. The zero-th order terms and include the protein synthesis and the release from the mitochondria processes, the linear terms include the protein degradation, chemical autoprocessing and protein complex dissociation; the quadratic terms take into account the activation and de-activation of protein P _i , the polinomial terms describe the protein association into larger complexes. No mass conservation constraint has been imposed to the system.

In our approximation we considered both the topology of the protein interaction map and the kinetic parameters as constant in time, i.e. each protein keeps the same neighbours during the time evolution of the system and interacts with them with constant strength. We decided to completely assign the connectivity matrix of the network on the basis of the existing experimental data. On the other hand, the kinetic parameters were largely unknown on the basis on the same information sources: as a consequence, in this application, the object of the "inverse problem" are the unknown model's constants. The "inverse problem" has been implemented with the following scheme:

1. eqs.(5–6) are solved and the time course of variables P _i and x _i (i = l...N) are calculated for a given set of model's parameters

2. the predicted time course of certain quantities is compared with the corresponding experimental data and a specific "distance" between time-courses evaluated

3. procedure is iterated up to minimizing that distance

Although, at least in principle, the strategy is simple, in practice the space of parameters to be estimated is very large, thus the strategy of points (1–3) above must rely on the availability of an efficient optimization algorithm. We have resorted to choose Genetic Algorithms (GA) for a number of reasons which will be highlighted in the following section.

The genetic algorithm (GA) is a programming technique that mimics biological evolution as a problem-solving strategy. Given a specific problem, the input to the GA is a set (called a "population") of potential solutions (called "individuals") to that problem. Each individual contains a "genome" able to provide a sub-obtimal solution to the problem. This ability could be quantified if a specific fitness function is defined, able to quantify how much an individual, by means of its genome, is fit for the solution of the optimization problem (i.e. to measure the "distance" between the sub-optimal and the optimal solution). The purpose of the GA is to produce successive population of individuals which are generated with the aim of increasing, as much as possible, the fitness of their individuals, i.e. their ability to solve the optimization problem by decreasing that "distance". This is done by producing successive populations of individuals by using the same procedures of the natural selection: mating and mutation. In the GA workflow, given an initial population of individuals, these are evaluated and classified according to their fitness. A selection rule is then defined to allow mating of couples of individuals, that mix their genomes, to form new ones (a further population) and an appropriate frequency of mutation of the genomes is defined, to introduce "new tracts" into individuals (which, in turn, would have been composed only by tracts coming from previous populations). If selection rules for mating and frequency of mutation are appropriatly chosen, the GA produces successive sets of individuals (" generations") which are progressively more and more fit to the optimization problem. In other words, individuals are better and better approximation of the optimal problem's solution.

The " inverse problem" we have attempted to solve starts from the description of a signalling network in terms of biochemical interacting species and reaction's constants. After a mining procedure to discover the value of the known reaction's constants, the system of eqs. (5–6) can be solved, by setting, for the unknown reaction's constants, an initial gauge of values. The solution of eqs.(5–6), in terms of functions describing the predicted time course of each of the system's variables (i.e. the concentration of all the biochemical species of the network), is thus strictly related to the intial set of reaction's constants. If one defines, as individual of the GA, the complete set of reaction's constants (the l known constants and the N - l unknown constants), its ability to produce an optimal solution to the problem can be measured by evaluating the " distance" between the predicted time-course (f _pred ) of some variables and that effectively measured by an experimental test of the same variables on that network (f _exp ). Formally, a distance between the two functions representing the j-th variable can be defined as follows:

where t is the (discrete) time length of the trajectories spanned by the variables. If one has k experimentally measured variables, the overall distance between that solution and the "optimal" solution would be

Eq.(8) can be thus retained as the "fitness" function of the considered individual; one can thus measure its "distance" from the "optimal" solution. Indeed, a more general formulation of the fitness function could be given by attributing "empirical" weight factors α to each variable, as to produce a different impact on the overall d value

The aim of the GA is to produce solutions which progressively reduce the value of the distance of its individuals. The scheme of producing successive "generations" of individuals can be resumed as follows:

1. start with a set of initial N individuals {K _i , i = 1, n}, each consisting of the same l known constants and by a number n - l of randomly selected guesses of the unknown constants (Fig. 2). Each value of K _i is a real number in the interval [10^-5,10⁰]. The interval was chosen on the basis of a reasonable number of kinetic values of protein-protein interactions published in the literature

2. for each individual, evaluate the distance d of eq.(8)

3. select, according to some defined rule, the individuals to be mated to form the new generation of individuals.

4. perform the mating procedure as follows: given two different individuals {K ^A (i)} and {K ^B (i)}, we randomly select the index m (l <m) and join the two individuals to produce a new individual {K ^A (i + 1)} such as

The parameter estimate does not include the topology of the network, that is the connectivity matrix is considered as a constant of the system and no interaction parameter is allowed to go to zero during the optimization procedure. The experimental data used as model constraints to optimize the system are the experimental time course of the concentrations of the active fraction of ERK-1, c-Raf, MEK, PKC-iota proteins [96–98]. These data were obtained measuring the phosphorylation level of these proteins by optical methods, following delivery of NGF to the cell. They are used to calculate the Fitness Function, upon which the GA is based. The optical signal were sampled every 2 minutes for 1 hour, that is a total of 30 sampling points for each protein were used to fit the system. The whole set of model parameters includes 278 indipendent values, of which 15 were extracted from the literature and the other were fitted using the GA.

The algorithm starts assigning every individual with random genome. The initial genes g _i were randomly generated in the following manner: g _i = 10^α where α comes from a flat (white noise) distribution in the range α ∈ [-5, 0]. This guarantees that the distribution of the initial parameters is flat in the logarithmic scale. The range is expanded in proportion to the number of components NC where the node (i) is a protein complex and if NC > 2: the range becomes α ∈ [-5 * (NC - 2), 0].

The Fitness Function F() is here defined, for each individual, as the inverse of the squared Euclidean distance between the experimental time course of the concentration of the activated fraction of ERK-1, c-Raf, MEK, PKC-iota proteins (see above) and the simulated time course for the same species, obtained using the genome {K₁,...,K _n } of the individual (Fig. 2, step 2) as parameters set; this distance is evaluated across the whole time interval (60 minutes), with a sampling time of 2 minutes:

Here p = p₁...p _np indexes the protein species used for the fitness evaluation, t = t₁...t _nt indexes the sampling times, (t) and (t) are the experimental and simulated time course of the concentration of the activated fraction for protein species p. The fittest genomes, those with the largest value of the Fitness Function, are given a greater probability to be selected to give birth to the next generation of individuals. The probability Pi of selection of the i ^th individual of the population is calculated as:

S _i = F^1/t, i ∈ {population}      (12)

where 10^-4 <T < 1 is a constant parameter, used to shape the distribution of the probabilities. It is worth underlining here that the best individuals tend to be selected at each generation, but the probability distribution gives any individual at least a small chance of being selected. The best individuals, then recombine their genes by the cross-over (Fig. 2, step 4), exchanging randomly selected but corresponding segments of the genomes, and eventually the offspring form the next generation (Fig. 2, step 5). The genes of the offspring are also allowed to randomly mutate with a low probability 0.005 <P _mut < 0.04. The individuals are distributed among NSp sub-populations, each containing NI of individuals, in our case NI = 16 and 7 < NSp < 33. The evolution process takes place independently within each sub-population at each generation. Every NM generations, with NM of the order of the sub-population size, MI of the best individuals in each population, again selected according to a probabilistic rule, move into a different sub-population, there replacing others that on their turn entered another sub-population: MI is of the order of 10%-30% of NI. This "migration" operator allows a sub-population to partially renew its genetic pool and tends to fasten the evolution process. The algorithm keeps in memory the "optimal" genome and the corresponding fitness, that is the best individual out of all the sub-populations obtained until that stage in the evolution process: these are compared with the best genome and corresponding fitness in the current generation: if the new fitness is better the optimal genome is replaced by the new one. The plot of the optimal fitness versus the generation number describes a monotone non-increasing function: when the curve derivative saturates, the procedure comes to an end and the individual corresponding to the optimal fitness provides the solution genome.

The GA in intrinsically parallel, thus the necessary computation can be very efficiently distributed over several CPUs. The GA was implemented on a cluster of Alpha CPUs, using the Fortran 90 language and the MPI protocol, under Linux operating system. In this implementation each computational node stores the genomes of a single sub-population, which evolves independently, except when there is a migration of individuals. In that case genome vectors are exchanged between the nodes (Fig. 2, step 4). In order to optimize the distribution of the computational load, data communications were reduced, which was exactly compatible with rarely occurring genome migrations.

In this work we have discussed the problem of mining, measuring and estimating the value of parameters needed in mathematical models describing the signalling processes mediated by protein-protein interactions. The lack of kinetic interaction rates measured in reliable in vivo and in vitro experiments is currently the major limitation to the creation of complex models of signaling pathway. We have attempted to show that biological information can be also extracted from a model which, levaraging on known kinetic parameters, attempts to provide a qualitative estimate of unknown parameters, even in the case of ill-conditioned optimization problems. We have thus sampled the space of model parameters using the Genetic Algorithm to estimate sets of unknown parameters. This sampling procedure has shown the existence of a basin of attraction for several kinetic constants. This might be interpreted as a a necessary condition for the network to produce a specific outcome of the time – course of its components. The estimated value of some of the parameters have shown a small coefficient of variation across the set of solutions, though the high dimensionality of this space allows to estimate realiable values and draw conlusions only on these few parameters.