WO2005013173A2 - Methode et systeme de selection de cibles therapeutiques par l'utilisation de reseaux dynamiques d'interactions moleculaires - Google Patents

Methode et systeme de selection de cibles therapeutiques par l'utilisation de reseaux dynamiques d'interactions moleculaires Download PDF

Info

Publication number
WO2005013173A2
WO2005013173A2 PCT/FR2004/002064 FR2004002064W WO2005013173A2 WO 2005013173 A2 WO2005013173 A2 WO 2005013173A2 FR 2004002064 W FR2004002064 W FR 2004002064W WO 2005013173 A2 WO2005013173 A2 WO 2005013173A2
Authority
WO
WIPO (PCT)
Prior art keywords
molecules
state
graph
vertices
network
Prior art date
Application number
PCT/FR2004/002064
Other languages
English (en)
Other versions
WO2005013173A3 (fr
Inventor
Todor Vujasinovic
Original Assignee
Helios Biosciences
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Helios Biosciences filed Critical Helios Biosciences
Priority to CA002534401A priority Critical patent/CA2534401A1/fr
Priority to EP04786022A priority patent/EP1649405A2/fr
Publication of WO2005013173A2 publication Critical patent/WO2005013173A2/fr
Publication of WO2005013173A3 publication Critical patent/WO2005013173A3/fr
Priority to US11/342,707 priority patent/US20060235670A1/en

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G16B5/30—Dynamic-time models
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G16B5/10—Boolean models
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G16B5/20—Probabilistic models

Definitions

  • the present invention relates to the field of integrative analysis of molecular interactions in a biological system. It relates in particular to methods for obtaining and analyzing networks of biological molecular interactions making it possible, from obtaining experimental data, to identify and describe the functioning of these interactions at the same time (i) between molecules interacting two by two, (ii) at the level of the results of the interactions exerted on a given molecule, and (iii) at the level of the whole network of interactions considered. Even more particularly, this method of analysis makes it possible, once the functioning of these interactions has been described, to predict the consequences, over the entire network of molecular interactions considered, of activations or inhibitions of the molecules forming this network. It thus makes it possible in particular to identify potential therapeutic targets, to understand the mechanisms of action of xenobiotics.
  • the approach which is increasingly used to develop new pharmacological treatments is another approach, called “physiological” or “comprehensive”, which consists in exploring and understanding the pathophysiological mechanisms, and in particular the mechanisms molecular pathophysiology, leading to the pathology concerned, in order to define the molecules of the organism to be treated, which will be the target molecules (or “therapeutic targets") of chemical treatments.
  • the identification of these target molecules then makes it possible, in a second step, to carry out sieves of potential therapeutic molecules of synthesis in order to identify those which will directly modify the biological activity of these therapeutic targets, or even to perform oriented syntheses of such therapeutic molecules when the spatial structure of the target molecules is known.
  • Integrative biology methods aim to analyze the role of molecules present in the organism to be treated, taking into account (and therefore integrating) in this analysis the other molecules with which they interact. Their objective is therefore to make it possible to obtain models of the cascades (or networks) of molecular interactions of living organisms, in particular those involved in pathological processes. In the context of the selection of therapeutic targets, such models aim to be applied to select these targets. More precisely, such an application must make it possible to predict the consequences of the actions of activation or inhibition of the molecules of the network, in order to identify those which will have a therapeutic effect. It is only conceivable to make such predictions on a large scale and in a sufficiently reliable manner if the model allows systematic simulations of the effects of inhibition or activation actions of the molecules of the cascade.
  • the modeling methods proposed today are, on the one hand, methods producing static models and, on the other hand, methods producing dynamic models.
  • the modeling methods producing static models consist in constructing static graphs representing cascades of interactions of biological molecules from data in the scientific literature (publications in journals, analysis of expression profiles of molecules, taking into account sequence data, etc.).
  • the resulting graph can be represented in the form of a diagram, most often in two dimensions, whose nodes (or vertices) of the graph are the molecules, and where these nodes are connected by lines or arrows (or arcs, or vertices of the graph) representing the interactions between molecules.
  • Examples of static graphs are those constructed in various public databases such as for example the KEGG database (M. Kanehisa and S. Goto: KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Research, 28 (1): 27- 30, 2000).
  • Diffusion the diffusion of molecules in the biological system studied or outside the biological system studied (for example a cell) is also taken into account, as equivalent to a synthesis or a degradation (respectively) within the system.
  • Variables used all the existing methods define the variables as being the rate, or the concentration, or the total quantity, of the molecules, noted here Xj for the molecule i, and not its proportion of variation compared to a standard state x i0 .
  • Continuous functions for a continuous function, the variables change continuously (as is the case in real biological systems) and not discrete.
  • Deterministic model once the model has been calculated, the network can only pass from one state to another by a single path (unique sequence of intermediate states). The fact that a model is deterministic makes it possible to obtain a linear growth in the quantity of calculations during simulations. Conversely, in non-deterministic models, the amount of computation required during simulations tends to increase exponentially with the size of the network, which can lead to an impossibility of implementation for large networks.
  • Level A knowledge of the existence in itself of molecular interactions, and at least part of their orientations and part of the effects of interactions (activation / inhibition or synthesis / degradation). Only level A knowledge is widely available to date. Therefore, only a method requiring only level A knowledge can be applied to wide area networks.
  • Level B level A with all the directions of the interactions and all the effects of the interactions.
  • Level C extensive functional knowledge of the network, i.e. level B plus other data such as: chemical reaction rate constants, description of threshold effects, description of allosteric effects, etc. To date, whatever the living organism considered, for most of the molecules in the molecular interaction networks, level C knowledge is not available. A detailed functional description of the biological network is necessary for the implementation of the method when level C knowledge is required.
  • any method requiring this type of knowledge for its implementation can only be applied to very small well-studied and known networks (a few dozen molecules at maximum) and is in fact unsuitable for its application to large networks (more than 100 to 150 molecules).
  • Level A linear growth with the size of the network (in number of molecules) of the required amount of computation. This corresponds to the possibility of implementation on a standard power server (general public). The methods implementing calculations, the quantity of which increases linearly with the size of the network, can be applied to large networks (provided that there is no other limit to this application). Level B: growth in the amount of intermediate computation between cases A and C. The methods implementing computations whose quantity increases in an intermediate manner between A and C are theoretically applicable to large networks but at a high or very high cost (and subject to not presenting any other limit to this application).
  • Level C exponential growth with the size of the network (in number of molecules) of the amount of calculation required. Any method implementing calculations, the quantity of which increases exponentially with the size of the network, requires very high computing power. For example, some applications of Bayesian networks require around 30 minutes of computing time on a server equipped with a 1.2 Giga Hertz processor for a network of 8 molecules: on a network of 32 molecules, the time to calculation on the same computer would in this case be more than a year and a half. In practice, even with today's most powerful computers, the methods exhibiting an exponential growth in computing time are not applicable to large or very large networks (a few thousand to a few tens of thousands of molecules and more; some of them are not applicable even to networks of a few hundred molecules). (3) Maximum size of network implemented: this is the maximum size of the networks on which the method has been successfully implemented to date.
  • the present invention aims to provide a method for obtaining dynamic models of molecular interaction networks in a biological system, which make this type of application possible. For the sake of clarity of this text, a certain number of terms are defined below.
  • This increase (or decrease, respectively) in biological activity can correspond either to an increase (or a decrease, respectively) in the number of molecules of a given type present in the biological system analyzed, each keeping the same activity (or biological function, either to an increase (or a decrease, respectively) in the activity of molecules of a given type, their number remaining constant, or to a combination of these two mechanisms, or to the result of these two mechanisms.
  • Activation (or inhibition, respectively) can also be: the consequence of an increase (or a decrease, respectively) in the number of molecules associated with a decrease - (or an increase, respectively) in their biological activity , if the overall result is an overall increase (or an overall decrease, respectively) in activity, and vice versa.
  • Activation can be non-zero or zero depending on the molecules considered and the biological system considered. It can be variable over time. The fact that certain interactions of the network of molecular interactions considered correspond to zero activation (or inhibition, respectively) is only a particular case of the field of the invention.
  • the biological activity of a biological molecule (s) considered corresponds to any capacity of the molecule (s) considered to have a chemical and / or physical interaction with any other molecule. another type (or with another molecule of the same type). This chemical and / or physical interaction may or may not result in the acquisition (or loss) by one of the interacting molecules of capacities to have a chemical and / or physical interaction with any other molecule of another type (or with another molecule of the same type).
  • Chemical interactions are any interaction between two (or more) molecules causing a chemical reaction (which can be represented by a modification of the formula of a molecule, or the synthesis, or the degradation of a molecule).
  • Physical interactions are any interaction between two (or more) molecules causing the formation of a stable or unstable complex between these molecules.
  • Examples of biological activities of molecules and corresponding molecular interactions are (not exclusively): the activity of activating the transcription of a given gene (molecular interaction: protein (transcription factor) - DNA), the activity of implementing a chemical reaction (molecular interaction: protein (enzyme) - molecule (substrate), allowing the transformation of the molecule-substrate into molecule-product of the chemical reaction), the activity of formation of '' a protein molecular complex itself having such or such biological activity (molecular interaction: protein (complex subunit) - protein (complex subunit)), etc.
  • biological molecule it is understood here any molecule, whatever its complexity, present in the biological system considered.
  • biological system it is understood here any living organism, whether it is prokaryotic or eukaryotic, and whether it is unicellular or pluricellular, and that the biological system corresponds to this organism as a whole or to a part of this organism.
  • biological system it is understood here any living organism, whether it is prokaryotic or eukaryotic, and whether it is unicellular or pluricellular, and that the biological system corresponds to this organism as a whole or to a part of this organism.
  • Whole organisms - A cell (eukaryotic or prokaryotic) as a whole. - A set of cells interacting directly or indirectly with one another, or not interacting with each other: 0 all of the cells in culture in a petri dish; 0 all of the cells forming an organ or part of that organ: the tonsil nucleus of a mammalian brain. - A multicellular living being. - The different examples plus their environment. Part of an organism: - An organelle of a cell, such as a mitochondrion.
  • a set of molecules participating in a given biological function such as a set of molecules participating in cellular respiration, or a set of molecules participating in cell death, whether this set of molecules is made up of all the molecules participating in said biological function where only a part of them.
  • the set of molecules forming the network of molecular interactions as described in the form of a static graph in FIG. 2 is an example of a biological system.
  • Many static graphs are, for example, available in the public KEGG database (M. Kanehisa and S. Goto: KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Research, 28 (1): 27-30, 2000). Every biological system is made up of molecules, these molecules interacting with each other in a more or less stable and variable manner over time and the effects of the environment of this system on the biological system itself.
  • apoptosis is the result of the interaction of multiple molecules (hormones, proteins, second messengers, etc.) which, for some of them, have physical interactions or chemical more or less stable over time.
  • network of molecular interactions is understood here the set of molecules analyzed by the method of the invention associated with the set (or part of this set) of their possible biological interactions.
  • the network can include all the molecules of the biological system concerned, or only a part of these molecules.
  • the network can be represented visually in the form of a graph (as an example is given in the description below). It is this type of visual representation that is behind the use of the term "network". However, such a representation is not a prerequisite of the invention.
  • the network can also be represented by a table (or a matrix) for which each row corresponds to one of the molecules of the network and whose columns correspond to the characteristics of the possible biological interactions of these molecules (or of a part of these interactions or their characteristics).
  • a graph is here a representation of the network of molecular interactions in the form of a graph whose vertices (or nodes) correspond to the molecules of the network of molecular interactions represented and whose edges (or arcs) connecting the vertices correspond to the interactions molecules of the network of molecular interactions represented.
  • vertices or nodes
  • edges or arcs
  • variable associated with a vertex of the graph it is understood here a quantitative variable in the mathematical sense of the term, which can take numerical values, and whose value at a given state of the graph represents the state of the corresponding vertex with regard to a quantity relating to a molecule of the biological system considered.
  • this quantity can be a level of expression of a gene expressed in the biological system (for example, the abundance of messenger RNA, measurable in particular by the microarray technique), a level of abundance of a protein, a level of activity of a protein, a level of abundance of a metabolite, etc., provided that the quantity considered is measurable experimentally, by a direct means or not.
  • a state of a graph is a graph for which a numerical value is given for each variable (associated with each vertex). The case where a nonzero numerical value is given only for a part of the variables (and associated with the corresponding vertices), another part of the variables (associated with other vertices) being zero, is only a particular case state of the graph.
  • a state of the given graph is a representation of a real or simulated state of the corresponding network of molecular interactions, and by extension a representation of a real or simulated state of the corresponding biological system.
  • giving a variable associated with a vertex of the graph the value zero can correspond to a representation of the situation where the molecule corresponding to this vertex is not present in the network of interactions (which does not mean that it is not present in the biological system), or else in the situation where its biological activity is zero.
  • the fact of giving a null value to a certain number of variables therefore corresponds to considering that at a given time these do not interact with the rest of the network, but their value can become non-zero at another time following a modification of the network status. Giving a variable a null value therefore does not necessarily amount to excluding the corresponding vertex from the network.
  • the invention also relates to a computer system for obtaining a dynamic model of a network of molecular interactions in a biological system, and the analysis of these molecular interactions when a stimulus is applied to the dynamic model, comprising at least at least one central data processing unit connected to at least one quantitative experimental database, the computer system comprising: A) a module for constructing a static graph, the vertices of which represent biological molecules and the arcs of which represent physico-chemical interactions existing between these molecules, each vertex being associated with a quantitative variable measured experimentally and each arc of the graph being associated with a mathematical relationship; and B) a learning module for calculating the parameters of each relation from quantitative experimental data concerning the vertices of the graph, by implementing learning techniques using gradient descent used for the configuration of networks.
  • the computer system according to the invention can also comprise: C) a simulation module for performing several iterative simulation procedures consisting in imposing a stimulus on a graph state measured experimentally and chosen as "state to modify", the stimulus modifying the value of one or more of the quantitative variables associated with the vertices of the graph, thus constituting a starting state of the simulation from which a propagation calculation is performed within the graph, to obtain a "final state of the graph ”; and D) an iteration module for modifying the stimulus.
  • the computer system according to the invention may further comprise: E) a module for calculating proximity between the “final state of a graph” and the “state to be modified”, or between the “final state of a graph ”and a desired state, and hierarchy of vertices and stimuli imposed on the vertices of the graph, the hierarchical vertices corresponding to classified therapeutic targets.
  • the computer system forms a tool for analyzing biological experimental data, and in particular a tool for hierarchy of biological molecules vis-à-vis a biological problem.
  • a first aspect of the invention is a method for obtaining a dynamic model of a network of molecular interactions.
  • the real sign of the term (ii) is determined by the result of the calculation of its parameter (s). This term (ii) is of opposite sign to term (i) once the parameters have been calculated, but this does not necessarily appear in its mathematical formulation, where the sign of the parameter (s) is not specified a priori shareholder (s).
  • each quantitative variable associated with a vertex represents the relative variation of the quantity of the molecule corresponding to said vertex, compared to a standard state of the biological system.
  • the "quantity of the molecule associated with a vertex" can relate to any aspect that can be measured directly or indirectly, this molecule, whether it is its concentration, its activity, its level of expression , etc.
  • said standard state is preferably a stable state of the biological system, in which the quantity of each molecule associated with a vertex of the graph is measurable experimentally.
  • this standard state can correspond to a given physiological state (for example healthy or sick) actually observable, or to an artificial state of the system, for example in the state a pool of several biological samples taken under different experimental conditions.
  • the relative variations in quantity of the molecules of the network are therefore represented in the form of variables dependent on the relative variations in quantity of the molecules interacting on them (ie in interaction with them and upstream in the network in terms of propagation of activations / inhibitions).
  • the definition of the variables corresponds directly to the experimental measurements available: in fact, in most molecular biology technologies (including screenings of messenger RNA expression), the absolute quantity of molecules present in the biological system of interest .- is not measured (nor measurable); only the proportion of their variation compared to a reference state is measurable.
  • the terms inertial (i) and return to the initial state (ii) make it possible to calculate the values of Xi and the variations of values of Xi over time as a function of the values of Xj (1- »n) and of the variations of values of X j (1-» n) over time.
  • inertial term is meant: - resistance to change, resistance in particular initial, and - a time to arrive at the maximum variation, which makes it possible to account for the complexities of propagations in the network.
  • the purpose of the inertial term (i) is to make it possible to integrate a resistance of the variables to change and a time difference between the modifications of the variables upstream and downstream of the network. It introduces in particular: - the integration of the time factor - the taking into account of the propagation speed differences within the network as a function of the sub-circuits - the taking into account of the temporal delays consecutive to the influences of the feedback loops on the propagation in the network, and - it makes it possible to calculate the kinetics of molecular interactions within the network directly from experimental data, without prior knowledge of the speed constants of these kinetics, and without making a priori on any other parameters.
  • This inertial term (i) tends towards a finite limit, which makes it possible to avoid significant divergences during simulations (improvement of their reliability): this avoids the risk of divergence (or "explosion") of the values of the linked variables to iterative propagations in feedback loops or during simulations involving extended times.
  • divergence or "explosion"
  • the formulation of this inertial term is preferably not very restrictive, and makes it possible to account for multiple forms of relations.
  • any initial variation in the effect of term (ii) on the calculated value of Xj can therefore be considered as consecutive to a prior variation in the effect of term (i) on the calculated value of Xj.
  • Xj can only present a variation if, at at least a given time, the variation of Xj at the following time calculated by the term (ii) is less in absolute value than the variation of Xj at the following time calculated by the term (i).
  • Xj can only present a variation, from a stable state, if, over at least a given time space, the variation of the calculated value of the term (ii) is less in absolute value than the change in the calculated value of term (i).
  • This characteristic is inherent in the fact that the term (i) is expressed as a function of X j (1 ⁇ n) while the term (ii) is expressed as a function of Xj.
  • the variation of term (ii) is initially lower in absolute value than the variation of term (i).
  • the effect of term (ii) on the variation of Xj may or may not become greater than the effect of term (i) on the variation of Xj. If this is the case, Xj will tend to return to its initial value.
  • the methods for obtaining a dynamic model of a network of molecular interactions in a biological system comprise a second step (step B), in which the parameters of the associated relationships are calculated. to each of the arcs of the graph, from quantitative experimental data concerning the vertices of the graph. This calculation is preferably carried out by the use of learning techniques.
  • a dynamic graph entirely deterministic, consisting of the static graph at the edges of which are now associated mathematical relationships whose parameters have all been defined numerically.
  • This calculation step can be performed by the use of learning procedures used for the configuration of artificial intelligence networks, for example those developed in computer science in "neural network” methods (including recurrent neural networks) by "simple” gradient descent (taking as a basis for calculation the pairs of data (Xj, X j ) provided by the experimental data independently of each other), or by gradient descent in time (where these couples are not considered independent).
  • the data pairs (Xj, X j ) provided by the experimental data are defined as follows: either i a molecule of the network, represented by the commit i, and either j any molecule of the network interacting on i, represented by the vertex j.
  • Xj and X j are the variables associated with vertices i and j, respectively.
  • experimental measurements of the values of Xj and X j under given experimental conditions and at given experimental times make it possible to obtain numerical values of Xj and X j .
  • a pair of experimental data (Xj, X j ) corresponds to the measured values of Xj and X j at a given experimental state (same time, same experimental condition).
  • the experimental data used to carry out step B) mentioned above have the following characteristics: Nature of the experimental data.
  • Those data are quantitative data concerning the molecules (corresponding to the vertices of the graph) and are for example levels of expression of genes expressed in the biological system (by measuring the abundance of messenger RNA, for example by the microarray technique DNA levels) and / or protein abundance levels and / or protein activity levels and / or metabolite abundance levels.
  • these data can be expressed in the form of a proportion of variation in quantity compared to a reference situation (standard state).
  • indexing of experimental data can be automatically indexed to the graph.
  • the role of this indexing is to link each experimental data to the corresponding biological object of the graph (vertex of the graph, or stop of the graph for the data pairs (Xj, X j ), so that these two types of information (experimental data and graph) can be used together when implementing the parameter calculation system.
  • Many commercial or free database systems allow this indexing to be created without any particular technical difficulty for those skilled in the art of biology or bioinformatics (Oracle commercial databases, Microsoft SQL server, FileMaker , open access databases: postgreSQL).
  • this indexing step may not be necessary per se.
  • the experimental data for the values of the pairs (Xj, X j ) are in the form of expression kinetics.
  • expression kinetics is understood here a set of series of experimental data ordered in time, each series of data corresponding to a set of values of couples (Xj, X j ) measured experimentally at a given time.
  • Each series of data can concern either the set of vertices of the graph, or only a subset of these vertices.
  • the different times correspond to successive times during the observation of a biological process using the biological system modeled by the graph, whether this process is natural or artificially induced in the laboratory.
  • Such kinetics preferably comprises at least three successive times, and, to improve the quality of the calculation of the parameters, more than three times.
  • step B) of the methods of the invention preferably takes account of the following principles:
  • the graph is considered to be in a stable state of reference at a given time, this stable state being measurable experimentally.
  • the stable reference state in question corresponds to an existing and measurable state of the biological system studied, which can be considered to be stable over time with respect to the modeled biological process.
  • the standard state which is arbitrarily defined by the experimental biologist, is used to carry out quantitative experimental measurements.
  • the stable reference state corresponds to a real state of the modeled system (Le., Not artificial), and serves as a reference for the calculation of the parameters of the model. It is considered a state of the system where the activation and inhibition processes at within the network are balanced, or have weak oscillations around a theoretical state of equilibrium. It represents the state towards which the system generally tends to return during simulations. It can be the same, or different, from the standard state.
  • the stable reference state is directly measurable experimentally as soon as the biological problem studied makes it possible to define a reference state of the biological system.
  • the measurement of the Xj of all the vertices of the graph in this state is used, in the calculation of the parameters, as a stable reference of the graph, in particular for the error minimization procedure.
  • the stable state is mathematically defined by the vector of the set of experimental values of the variables of each vertex measured in the corresponding biological state (measurements made for all the vertices of the graph).
  • the standard state for the measurements is the stable state.
  • the Xj are close to 1 in general (if the standard measurement state is time t 0 ).
  • V i, X, 1, and to introduce (in the sense of "adding") this vector of Xj at the initial kinetic time without having measured it.
  • the standard state is considered stable, arbitrarily. This is often possible if the standard state does not correspond to a pool of different biological tissues.
  • the experimental data are measured during kinetics (see above). In the case where the biological process of interest is studied during the transition from a state initial stable to a final stable state, and where experimental measurements are made at these two states and at intermediate times, two stable states are defined: the initial state and the final state of the kinetics of the biological process studied. However, having experimental measurements corresponding to two stable states is not a prerequisite for implementing the invention.
  • the calculation of the parameters of the relationships (Xj, X j ) is carried out by the implementation of learning techniques used for the configuration of artificial intelligence networks (such as those implemented for neural networks), from quantitative experimental data concerning the variables of the graph.
  • this calculation can use algorithms for digital resolution of propagation or back propagation with calculation of the error.
  • a second aspect of the present invention relates to a method for analyzing a network of molecular interactions in a biological system, comprising the following steps: A ′) use of a dynamic model of the network of molecular interactions, said model being likely to be obtained by a process described above, and constructed from a static graph whose vertices represent biological molecules of the biological system and the edges represent physicochemical interactions between these molecules, and from experimental data regarding the levels or activities of these biological molecules.
  • the propagation calculation within the graph can be performed for a number of time steps such that the duration of the simulation does not exceed the duration of the biological process to be simulated defined in step C). However, it is also possible to let the simulation continue beyond the duration of the biological process to be simulated defined in step C), for example if we are looking to see if the network will eventually find a new stable state (equilibrium state) and if we do not know a priori how long it will take. It is important to note that the duration of the simulation defined in step C) can be longer than that of the experimental kinetics used for the calculation of the parameters (or shorter).
  • steps C), D) a) and D) b) above are carried out, using (without reconstructing) a dynamic model of the network of molecular interactions chosen, said model being capable of being obtained by a method such as the methods for obtaining dynamic models of networks of molecular interactions described above.
  • Another particularly important aspect of the present invention is a method for selecting therapeutic targets implementing a dynamic model of a network of molecular interactions in a biological system, by implementing a computer system, comprising the steps and following characteristics:
  • D) several iterative simulation procedures are carried out, each comprising the following steps: a) a stimulus is imposed on the state to be modified, that is
  • step D) b) can be continued beyond the duration specified in step C)
  • Step A ') of the above methods can be carried out in the same way as steps A) and B) of the methods for obtaining dynamic models of interaction networks described above.
  • step C) can be carried out by taking account of the following elements, when the case permits:
  • Time step The duration of the biological process to be simulated is divided into a series of time steps, regularly spaced or not; the time steps are defined so as to be preferably smaller than the real experimental durations separating the series of quantitative experimental data used for the computation of the parameters of the relations.
  • the definition of these time steps is made necessary by the fact that any computer process of dynamic simulation consists in calculating states at discrete times, making the discretization of the time necessary. We therefore obtain a series of consecutive times, on which the simulation will be performed.
  • the first beat in the series is called the initial beat. This initial time corresponds to the starting state of the graph, defined below.
  • the state to be modified and the state to be achieved are defined as follows: We simulate actions on the vertices of the graph from a state to modify the graph identical or similar to its state as observed experimentally in the pathological condition (for example by screening the expression of messenger RNA on chips DNA from pathological tissue).
  • the state to be reached is defined as being a state close to a reference non-pathological state (as also measured by the experimental observation of the non-pathological condition, for example by screening for expression of the messenger RNAs on microchips). DNA from healthy tissue).
  • the simulation process then consists in identifying the vertices, and the stimuli on these vertices, which, starting from the state to be modified (the pathological state), best allow the graph to evolve (in part or entirely) towards a state close to the state to be reached (non-pathological state).
  • the state to be modified is defined as above, and the state to be reached is defined as the state, or a state close, to that obtained experimentally during the administration of this treatment (as measured for example by screening the expression of messenger RNAs on DNA chips from pathological tissues which have been subjected to the treatment concerned).
  • the simulation process then consists in identifying the vertices, and the stimuli on these vertices, which, starting from the state to be modified (the pathological state), best allow the graph to evolve (in part or entirely) towards a state close to the state to be reached (pathological state under treatment).
  • This particular implementation can also be carried out by defining the state to be modified as any possible state G of the biological system studied (for example the healthy state), and the state to be reached as the state obtained after administration of the treatment concerned with the biological system in state G.
  • step D) is carried out by considering the following elements:
  • Stimulus A stimulus is imposed on the state to be modified. This stimulus is exerted in the form of the variation of the value of one or more variable (s) of the graph (corresponding to one or more vertex (s)), that is to say an increase or d 'a decrease in this or these value (s), depending on the desired simulation. The values of all other variables remain unchanged. We therefore obtain a new state of the graph, which is
  • starting state of the simulation.
  • the starting state and the state to be modified therefore only differ by the value of the modified variable (s), all the values of all the other variables being identical
  • This state is defined as corresponding to the first step of the simulation
  • the stimuli relate each time to a single vertex.
  • Propagation From the starting state of the simulation, a propagation calculation is performed within the network. This propagation consists in calculating the new values of all the variables in the next time step, resulting in a new state of the graph, and in starting again the calculation starting from this new state for the next time step, and so on. This propagation is prolonged during the number of time steps (thus the biological duration) defined by the experimenter according to the biological question posed. It can possibly be extended until the appearance of a new stable state of the graph (a new state of equilibrium), or be stopped before. At the end of this simulation, a new state ("final state”) of the graph is obtained.
  • the previous process is repeated with a new stimulus, relating to one or more other vertex (es) of the graph, or possibly relating to the same vertex (s) of the graph with the imposition a new value to the variable (s).
  • step E) is a hierarchy of the vertices, and of the stimuli imposed on these vertices.
  • This classification therefore corresponds to the classification of the vertices, from that on which a stimulus is most likely to result in the desired state from the state to be modified, to that on which a stimulus is least likely to have this effect.
  • Each proximity corresponds in fact to one and only one vertex and one and only one stimulation value on this vertex. If the desired effect is the improvement of a pathological state, this classification is that of potential therapeutic targets, most likely to least likely.
  • step D) c) the proximity of each final state obtained in step D) b) can be calculated either with respect to the state to be modified chosen in step C), or with respect to another state, measured experimentally or determined arbitrarily, and defined as P "state to reach", which can be, for example, a healthy state. It can be the reference state defined above.
  • step E) consists in classifying all the bytes (vertex (s) of the graph - stimulus) in hierarchical order (ascending or descending) corresponding directly to the hierarchical order (ascending or descending, respectively) of the proximities associated with them.
  • the vertices of the graph directly correspond to the molecules of the biological network, which are therefore hierarchized in fact.
  • the hierarchical classification of the molecules of the biological network corresponds directly to the hierarchical classification of these molecules considered as therapeutic targets.
  • the invention therefore makes it possible to obtain a classification of potential therapeutic targets hierarchized according to objective statistical criteria, as a function of experimental data (measurements of Xj) and of functional knowledge of the network (existence of molecular interactions).
  • the best potential therapeutic targets are considered to be those corresponding to the best proximity to the state to be reached.
  • this invention makes it possible not only to identify therapeutic target molecules, but also to predict the direction of the therapeutic action which it will be necessary to apply to these molecules (activation or inhibition).
  • the therapeutic targets are therefore selected from data concerning all of the molecules studied, and not only those specifically concerning target molecules, since the criterion used for prioritization depends on the evolution of the graph as a whole, therefore on the set of experimental measures of expression and / or activation of all the molecules represented in the graph , and not the simple evolution of experimental measures of expression and / or activation of only target molecules. It is therefore indeed an integrative method responding to current needs as defined above, in particular as regards diseases with multi-factorial determinism, clearly bringing progress compared to the methods of selection of existing therapeutic targets.
  • the target identification method described above has the following advantageous characteristics: -
  • the calculations are based on a non-probabilistic method, which eliminates any limitation in terms of calculation time, unlike the methods of stochastic equations and Bayesian networks .
  • the invention integrates quantitative experimental data, which differentiates it from qualitative methods (Boolean networks, generalized logical formalisms, formalisms based on rules), makes it possible to avoid constraints and hypotheses on the functioning of the network, and makes it possible to increase the reliability of simulations.
  • the fact of defining the variables as similar to the actually measurable experimental data makes it possible to calculate the parameters of the relationships in an optimal manner (without having to extrapolate the values of the variables).
  • a first hierarchical classification of the vertices, considered individually, is obtained by performing steps A) to E) by imposing, for each of the simulations of step D), stimuli that relate to a single summit; a step D2) is then carried out, corresponding to step D) in which the stimuli imposed on each simulation are exercised on two vertices, either by testing all the combinations of two possible vertices, or by limiting these calculations to the combinations of two vertices among a certain number of the vertices best classified in step E).
  • a step E2) of hierarchical classification of the associations of two vertices on which stimuli are most likely to have the desired effect is carried out from the set of statistical proximities calculated in step D2).
  • steps D) and E) can be repeated iteratively, each time increasing the number of vertices on which the stimuli are exerted.
  • the method can include a step D3) following step E2) and corresponding to step D) in which the stimuli imposed on each simulation are exerted on three vertices, either by testing all the combinations of three possible vertices, or by limiting these calculations to the combinations of three vertices chosen from a number of the highest classified vertices in step E) and combinations of two best classified vertices in step E2), said step D3) being followed by a step E3) of hierarchical classification of the associations of three vertices on which stimuli are most likely to have the desired effect.
  • Steps D4) and E4), with stimuli on 4 vertices can then be added, and so on.
  • These methods of selecting therapeutic targets preferably include a final step of statistical classification of the proximities of graphs of all the simulations carried out, integrating all the classifications previously obtained.
  • the stimuli exerted on these different vertices can be applied simultaneously or not, that is to say that the simulation can take account of a time difference between the stimuli exerted on different vertices.
  • the relationship between Xj and Xj is established, for at least part of the physicochemical interactions between the molecules of the network, by an inertial relationship arising from that of the harmonic oscillator in physics , associated with a sufficiently large damping factor to limit the oscillation to a single cycle.
  • (dXj / dt) corresponds to the term of return to the initial state ( ii)
  • Xi is a variable associated with the molecule i
  • dXi / dt is the derivative of Xj as a function of time d 2 X
  • / dt 2 is the second derivative of Xj as a function of time
  • X j is a variable associated with the molecule j
  • ⁇ ii represents the inertia of i
  • ⁇ y governs the return to the equilibrium state of Xj
  • the pulsation ⁇ y corresponds to the response time of Xj to the variation of X j
  • Wj j is a coupling factor representing the strength of the interaction between molecules i and j, corresponding to a weighting of the effect of each molecule j on molecule i vis-à-vis the result of the set of combined effects of all molecules j exerting an eff and on i.
  • the relationship between the variables Xj and X j , two by two is established by a sigmoid relationship comprising a delay factor associated with a linear decay function.
  • the present invention also relates to a method for determining the mode of action of a xenobiotic, consisting in implementing a method for analyzing a network of molecular interactions in a biological system, such as those described above, under the following conditions: (i) the biological system in which a network of molecular interactions is studied is affected by the action of xenobiotic; (ii) the state to be modified "chosen in step C), corresponds to a state observed experimentally before the administration of said xenobiotic; (iii) the modifications to be made during step D) are identified for that the calculation carried out in step D) b) shows an evolution of the system towards a state close to the state observed after administration of the xenobiotic.
  • Another aspect of the invention is a method of predicting undesirable effects of treatments applying a dynamic model of a network of molecular interactions in a biological system, by implementing a computer system.
  • steps A to E we analyze on parts of the graph representative of known physiological functions, by simulations using the same method as in the previous aspects of the invention (steps A to E, possibly A to Ek when steps D and E are repeated iteratively by applying the stimuli on combinations of vertices up to k vertices), the consequences of the application of the treatment on these target molecules.
  • This analysis consists in identifying the possible evolutions of these sub-parts of graphs towards new states close to other reference pathological states (as defined by the experimental observation of these pathological conditions, according to methods similar to what is described. upper).
  • observation during simulations of the evolution of the apoptosis subgraph towards a final state having close proximity to a reference state of this graph corresponding to an activation of this physiological pathway makes it possible to predict an effect of cellular toxicity of the treatment in the tissue or tissues concerned.
  • This aspect of the invention therefore consists in implementing an analysis method as described above, under the following conditions: (i) the biological system in which a network of molecular interactions is studied is concerned by the treatment; (ii) the modifications of step D) a) correspond to the modifications of the levels or activity levels of the target molecules observed or desired during the application of the treatment; (iii) step D) b) of calculating the evolution of the biological system is followed by an analysis of sub-parts of the system corresponding to known physiological functions, in order to identify any evolutions of these sub-parts towards states close to reference pathological states.
  • the present invention also relates to a method for prioritizing potential therapeutic targets for a pathology, comprising identifying therapeutic targets by a method according to the invention, then predicting the possible undesirable effects of a treatment targeting these targets, and finally to determine the "therapeutic benefit / adverse effects" ratio of an action on each of the potential therapeutic targets.
  • the present invention in its various aspects, is that it makes it possible to work on graphs or networks of interacting molecules comprising a large number of molecules.
  • the number of variables Xi of the network of molecular interactions considered is therefore preferably greater than 100, greater than 200, or even greater than 300.
  • the invention also relates to an analysis method as described above using the molecular interaction networks of the invention, said networks being associated to form a hypergraph of networks.
  • the number of variables Xi of each network of molecular interactions is less than approximately 100 and the number of networks associated to form the hypergraph is between 2 and approximately 100.
  • Another aspect of the invention is a method of extending the graphs from results of experimental screening for variations in the rate of expression or activity of molecules, applying a dynamic model of a network of molecular interactions in a biological system, by implementing a computer system.
  • the steps and characteristics of the method are the same as previously, the only modification consisting in the following adaptation:
  • the method is used to identify new molecular interactions. This can be achieved by coupling the method of the invention described above, with statistical methods for finding correlation between points in an n-dimensional space (for example factor analysis, hierarchical classifications, etc.) such as ( but not exclusively) those used to date to search for correlations of gene expression from the results of screening for messenger RNA on DNA chips ("clustering" of genes).
  • cluster 3.0 developed by the Laboratory of DNA Information Analysis of Human Genome Center, http: //www.ims.u-tokyo.ac.ip/imswww/index- e.htmllnstitute of Médical Science, Universitv of Tokyo, au Japan (4-6-1 Shirokanedai, Minato-ku, Tokyo 108-8639 JAPAN).
  • the "cluster 3.0" software is available on the website http: //bonsai.ims.u- tokyo.ac.jp/ ⁇ mdehoon/software/cluster/.
  • the experimental data used can, for example, be that produced by messenger RNA expression screens on DNA chips. This coupling consists in using the parameters calculated by the implementation of the invention to recalculate a new matrix of experimental data for measuring the expression of the rate or activity of the molecules, by eliminating matrices of experimental results from originate the molecular interaction factors included in the parameterized dynamic model (such as the dynamic or inertial resistance component), then carry out correlation searches.
  • This "cleaning" of the original result matrices consists in other words in eliminating the "statistical noise” linked to these factors, these factors then being considered as introducing distortions, in the measures actually observed of the expression rates. or activity of molecules, compared to what these measures would have been, from a theoretical point of view, in the absence of these factors.
  • the dynamic resistance of the expression of a given gene A (therefore the inertia of the modification of the level of corresponding messenger RNA) to two distinct stimulations exerted by the molecules B and C (themselves distinct ) can vary, preventing before any "cleaning" of this type from highlighting both a correlation between the expression of A and the activity of the molecule B, and a correlation between the expression of A and the activity of molecule C.
  • the invention therefore relates to the use of a dynamic model of a network of molecular interactions in a biological system capable of being obtained by a method as described above, to extend a graph. static whose vertices represent biological molecules and the arcs represent physico-chemical interactions between these molecules.
  • FIG. 1 represents a synthetic diagram of the various steps of a method for identifying targets according to the present invention.
  • the text within each rectangle corresponds to the abbreviation of the name of the protein as described in the text.
  • the parameters of the relationships between the vertices of the graph and the simulations were calculated on the whole of this graph (example 4).
  • Figure 3 shows graphical examples of calculated (triangles) and observed (squares) kinetics for some genes (example 4).
  • Figure 3A ORF YBL015W (ACM).
  • Figure 3B ORF YMR169C (ALD3).
  • Figure 3C ORF YIL125W (KGD1).
  • Figure 3D ORF YNL071W (PDA2).
  • Figure 3E ORF YAL054C (ACS1).
  • Figure 3F ORF YFL018C (LPD1).
  • FIG. 4 represents a diagram of the graph constructed in example 5. This graph includes 133 molecules in the network of molecular interactions (133 vertices), and 407 molecular interactions between these molecules.
  • Figure 5 shows examples of parameterization curves, in which the kinetics measured experimentally are represented in white and the kinetics calculated by simulation are represented in black, for some molecules (example 5).
  • Figure 5A ICL 1 (YER065C).
  • Figure 5B IDH1 (YNL037C).
  • Figure 5C ACM (YBL015W).
  • Figure 5D PCK1 (YKR097W).
  • FIG. 6 represents a classification scheme of the molecules of the network by hierarchical classification of the distances calculated between on the one hand the state to be reached and on the other hand the states obtained by simulation
  • the ordinates correspond to the calculated distance values.
  • the 133 molecules of the network are classified from left to right from that associated with the shortest distance to that associated with the highest distance, each point corresponding to a molecule in the network.
  • i be a given molecule of the network, and Xj its quantity (or its concentration) within the biological system studied.
  • x i0 be the experimental measurement actually carried out from i to a "standard state" of the system biological, used during measurements.
  • XH be the experimental measure actually carried out from i at an instant t. The variable used is:
  • the standard efaf is a measurable condition used to perform biological measurements, against which all other measurements are quantified. It can correspond to an artificial state of the system, for example to a pool of several biological samples taken under different experimental conditions (artificial state), or to a really observable state (not artificial) of the system. This variable corresponds well to the type of biological measurements actually achievable.
  • the measurement actually carried out for each RNA at a given experimental time t is the ratio of the signal emitted by the hybridization from RNA present in the biological sample at time t on the signal emitted by RNAs of the same type present in the sample at a standard state of the biological system studied (for example the initial time of the biological experiment).
  • RNA molecules not being directly measurable because it depends on experimental parameters not directly controlled (yield of the labeling reactions of the probes, yield of hybridizations on the chips, etc. , these parameters differing in a non-predictable way between two RNAs of different types given).
  • the quantity of signal measured in the standard state therefore serves as a measurement standard for that at other times, based on the assumption that for a given type of RNA, the experimental parameters influencing the signal finally emitted are the same.
  • Xi therefore corresponds directly to the quantitative biological measures that are really productive in the current state of molecular biology techniques.
  • Variables Xj, Xj etc. are therefore equal to (xit / xio), (jt / jo) etc. - (2)
  • n be vertices ji, j 2 , ..., jn of the graph (corresponding to n molecules of the network) which act on a vertex i (orientation of the graph from j to i).
  • Inertial term of these relations corresponds to a continuous function of the Xj. This term has an inertial component.
  • inertia we mean the fact that Xj presents a resistance to change following a variation of Xj: more precisely, this term of the relation must make it possible to account for the behavior according to variables: following a given variation of one or more of X j , the speed of variation of Xi will be initially low, then accelerate gradually.
  • the relationship between Xj and X j is established by an inertial relationship arising from that of the harmonic oscillator in physics associated with a damping factor large enough to limit the oscillation to a single cycle .
  • Xi variable associated with the molecule i dX
  • / dt derivative of Xj as a function of time d 2
  • Xj / dt 2 second derivative of Xj as a function of time
  • Xj variable associated with the molecule j
  • mj inertia of i (resistance to change)
  • ⁇ y damping (governs the return to the equilibrium state of Xj).
  • ⁇ pulsation (response time of Xj to the stimulus represented by X j )
  • wy coupling factor representing the strength of the interaction between molecules i and j, corresponding to a weighting of the effect of each molecule j on the molecule i with respect to the resultant of the set of combined effects of all the molecules j having an effect on i.
  • Xi; dXj / dt; d 2 Xj / dt 2 and Xj are data supplied experimentally or directly calculated from these data. For example, this data can be obtained by screening for mRNA expressions.
  • the relationship between Xj and the set of Xj is defined by a sum whose weighting includes a variable term over time and explicitly taking into account the respective velocities of the modifications of Xj, velocities represented by the derivatives: dX j (t) / dt, denoted in the sequence dX j / dt.
  • dX j (t) / dt velocities represented by the derivatives: dX j (t) / dt, denoted in the sequence dX j / dt.
  • ay is directly calculable, from experimental quantitative data, for each experimental time. This weighting factor varies over time.
  • the relationship between Xj and the set of X j is defined by a sum whose weighting includes a variable term over time and which takes explicit account of the respective accelerations of the modifications of the Xj, accelerations represented by the second derivatives: d 2 X j (t) / dt 2 denoted in the following d 2 X j (t) / dt 2 .
  • d 2 X j (t) / dt 2 denoted in the following d 2 X j (t) / dt 2 .
  • the relationship between Xi and the Xj is established by a sigmoid relationship comprising a delay factor associated with a linear decay relationship.
  • the relation associated with the edges corresponding to an inhibitory molecular interaction could be the opposite of that associated with the edges corresponding to an activating interaction (sigmoid curve decreasing or increasing, respectively), this characteristic being directly obtained during the calculation of the parameters. , by their positive or negative signs.
  • K ⁇ . [1 / (1 + e- ⁇ w0J - bl )] corresponds to the inertial term
  • Xi variable associated with vertex i.
  • Xj variable associated with vertex j.
  • wij weighting factor of the effect of j on i when several vertices (ji. J 2 , • Jn) act on the vertex i.
  • bi delay factor
  • the weighted result of the combination of the effects of the set of Xj on Xj is included from the start.
  • the combined effect of all X j is also represented here by a weighted sum, in the term inertial.
  • the relationship between the Xj and the X j is established if ilary to the previous aspect, but with the introduction of an additional weighting factor ay as defined by equations (4) or (7) above.
  • the parameter ay is a variable term over time and takes explicitly into account the respective velocities of the modifications of X j , or the respective accelerations of Xj, depending on whether a is defined by equation (4) or equation (7) , respectively. Since the times at which the experimental measurements are made are long compared to the accelerations, it is preferable define the weighting factors ay with respect to the velocities (equation
  • Equation (4)
  • Equation (7)
  • E ⁇ I [Xii (t) -X 2i (t)] 2 .dt it with: Xij (t): value of Xj at time t calculated by simulation with: X 2 j (t): value of X; at time t measured experimentally.
  • V square root of the term in square brackets: []. X ⁇ (t) and X 2 i (t) being defined as above.
  • the graph (or network) is entirely determined by the mathematical relations associated with the edges of the graph: it is a fully deterministic network.
  • the corresponding graph is oriented.
  • the network can be represented explicitly by the implementation of techniques for representing neural networks used in artificial intelligence. It is a non-Boolean network, neither Bayesian, nor organized in layers, allowing to represent redundancies of circuits and feedback loops. This deterministic network allows the implementation of simulations without significant computation cost, even for a very large graph.
  • propagation methods include the Neural Network Toolbox 4.0.2 software, developed in the computing environment
  • Propagation is inherent in the learning methods mentioned in Example 2.
  • the simulation step therefore consists in using the propagation method implemented in step B, or any other propagation method deemed adequate by humans. art.
  • thresholding procedures are associated with the propagation method chosen, in order to reduce the divergences (therefore to improve the convergences, that is to say the reliability). These can relate to:
  • a lower threshold i.e. any value of a predicted variable below 0 is reduced to 0 (or to a minimum background noise value if this can be defined by experimental data):
  • a higher threshold by imposing a maximum threshold on the values of Xj which can be obtained during simulations (for example by fixing a maximum threshold corresponding to a multiplicative factor (which can not exclusively be fixed between 1 and 10) of the maximum values observed experimentally for each Xj; this factor can possibly be defined during simulations of real experimental results available by testing several values of this factor -
  • the introduction of constraints in the loops during simulations This can be achieved by several methods, not exclusive some of others, this list is not exhaustive.
  • these simulations aim to analyze the biological system as a whole, and therefore to respond to the issues mentioned above.
  • These simulations therefore consist in describing the evolution of all the molecules in the network of molecular interactions over time following the initial "virtual” stimuli, including, if it turns out to be biologically important, until a new state balance of the graph.
  • These “virtual” stimuli can be applied systematically to each vertex of the graph, but it is also possible to carry out several stimuli at the same time, or sequentially, on several vertices.
  • Proximity of the states of the graph A statistical calculation of proximity between each final state calculated and the desired state, or between each final state and the state to be modified, is carried out. It allows, for each vertex, to associate a statistical criterion (proximity obtained) with the vertex on which the stimulus was exerted and with the stimulus exerted on this vertex.
  • the distance between two states of a graph can be calculated as follows:
  • a static graph corresponding to a network of molecular interactions in yeast was constructed by manual entry in a flat file txt type (without implementing a system automated), using data from the KEGG database: Kyoto Encyclopedia of Genes and Genomes, free access data at the internet address http.7 / www. ⁇ enome.ad.ip / keqq / keqq2.html.
  • This graph is more particularly centered on the mechanisms of cellular respiration (glycolysis, neoglucogenesis, metabolism of Pyruvate and acetyl CoA, etc.). It includes 116 molecules, enzymes or transcription factors. It includes 329 uni- or bi-directional interactions between these molecules.
  • each rectangle represents a protein.
  • the letters in the rectangle are the usual abbreviations of the name of the protein, according to the nomenclature of KEGG and SGD: Saccharomyces Genome Database developed by the Department of Genetics at the School of Medicine, Stanford University, USA.
  • This table presents the static graph data in a form which can be directly used by a computer modeling and simulation system as described in the invention.
  • the graph represents the 329 interactions between the 116 molecules in the network. Interactions are represented between the molecules two by two.
  • ORF codes are unique for a given protein and allow it to be identified without any ambiguity. They also make it possible to establish an unambiguous link with the results of screening on DNA chips (correspondence of the nucleic sequences of the corresponding messenger RNAs).
  • This publication describes an experiment in the culture of yeasts under conditions where the concentration of glucose in the culture medium gradually decreases (due to its use by yeasts for fermentation, since glucose is not added to the culture at any time. experience). Over time, the yeasts change their metabolism, their respiratory system going from functioning in fermentation to functioning in aerobic respiration.
  • This yeast culture has been studied over time, in particular by the practice of expression screening of almost all of the yeast messenger RNAs on DNA chips. These screenings were carried out at successive times, the results therefore producing an expression level kinetics for each messenger RNA.
  • the results show variations in the level of expression of a certain number of messenger RNAs over time, these being more particularly numerous among the messenger RNAs of cellular respiration proteins, a large part of which is represented in the graph described above. Under these experimental conditions, the graph that we have constructed therefore presents a dynamic evolution over time, which is represented by the kinetics of the molecules of the network (therefore the vertices of the graph).
  • each line corresponds to a molecule of the graph
  • the first column identifying the ORF (open reading frame) by its SGD code
  • Tables 6 to 8 below give examples of data for some molecules of the graph, data extracted from the page http://cmgm.stanford.edu/pbrown/explore/array.txt.
  • Tables 6 to 8 Examples of DNA microarray screening data for 2 of the molecules of the graph, from the results of the article: DeRisi JL, lyer VR, Brown PO: Exploring the metabolic and genetic control of gene expression on a genomic scale, Science. 1997 Oct 24; 278 (5338): 680-6. Full data are available on the web page: http://cmgm.stanford.edu/pbrown/explore/array.txt. Given the size of the complete data table, only part of it is shown here. Those skilled in the art will very easily be able to recover the data corresponding to the other molecules of the graph used in this example, on this internet page which corresponds directly to the table of all the data.
  • the expression screens of the messenger RNAs were carried out every two hours, for 12 hours, which corresponds to 7 experimental times (the initial time plus the 6 following times). These correspond to notations 1 to 7. The reader will find all the explanations corresponding to obtaining these measurements in the article cited in reference.
  • Name abbreviation of the name of the gene (according to SGD)
  • ORF open reading frame code (according to SGD)
  • G experimental condition corresponding to the "standard state" of the graph as described in the invention. This standard state, in this series of experiments, corresponds to the initial state of yeast culture. G1, G2, G3, G4, G5, G6,
  • G7 all correspond to the same standard biological sample.
  • R states of the graph at various experimental times.
  • R1 corresponds to the initial state of culture of the yeasts (same biological sample as G1), at time TO.
  • R2 TO + 2.5 hours
  • R3 TO + 4 hours
  • R4 TO + 6 hours
  • R5 TO + 7.5 hours
  • R6 TO + 9.5 hours
  • R7 TO +
  • G2-Bkg measurement of G2 minus the background noise during experimental measurements.
  • G3-Bkg measurement of G3 minus the background noise during experimental measurements.
  • G4-Bkg measure of G4 minus the background noise during experimental measurements.
  • G5-Bkg measurement of G5 minus the background noise during experimental measurements.
  • G6-Bkg measurement of G6 minus the background noise during experimental measurements.
  • G7-Bkg measurement of G7 minus the background noise during experimental measurements.
  • the variations in the measured values are linked to the variations in the yields of the various reactions used in the measurement method (DNA chips) and justify the use of a standard state as a measurement reference.
  • R1-Bkg measurement of R1 minus the background noise during the experimental measurements.
  • R2-Bkg measurement of R2 minus the background noise during experimental measurements.
  • R3-Bkg measurement of R3 minus the background noise during the experimental measurements.
  • R4-Bkg measurement of R4 minus the background noise during the experimental measurements.
  • R5-Bkg measurement of R5 minus the background noise during the experimental measurements.
  • R6-Bkg measurement of R6 minus the background noise during experimental measurements.
  • R7-Bkg measurement of R7 minus the background noise (experimental noise).
  • Steps A and B were implemented as follows:
  • the parameters were calculated from the experimental data by a classic method of learning neural networks, more precisely from the algorithms of the Runge Kutta method of back propagation through time (BPTT: back propagation through time ). The calculations were performed in double precision.
  • the effectiveness of the method was checked. Indeed, the average divergence result (global relative error) of the simulations compared to the experimental data is approximately 0.30, this divergence being essentially concentrated on 8 vertices (molecules) of the graph, for this series of data, whereas for the remaining 108 vertices (molecules), the divergence is very small.
  • This overall relative error result shows that the kinetics calculated during the simulations are close to the real data, since random kinetics would have given an error calculation greater than 1.
  • V square root of the term in square brackets [].
  • X-ii (t) value of X calculated at time t of the simulation
  • X 2 i (t) value of X measured experimentally at time t.
  • ⁇ sum of values at different times t
  • the rate of non-reproducibility of the experimental data can be estimated by the ratio R1. Ratio of experimental data (Table 8), and is overall 14% in this example.
  • table 9 gives the details of the calculations of the divergences of the set of kinetics during the simulations for all of the molecules of the network, in the form of a summary table.
  • Xii (t) value of Xi calculated at time t of the simulation
  • X 2 i (t) value of Xi measured experimentally at time t
  • This calculation amounts to calculating the difference in integral between the curves of the kinetics observed experimentally and the kinetics calculated during the simulations. It therefore concerns both the whole kinetics and the final state.
  • this variant does not of course consist in removing from the graph the "input" vertices, but in imposing that their kinetics remain the kinetics measured experimentally, this only during the simulations practiced during the error calculations of the procedure for calculating parameters by back propagation over time.
  • the parameters of the mathematical relations connecting these vertices to other vertices of the graph are therefore calculated, as for all the other edges of the graph, and the dynamic model finally obtained includes these vertices.
  • FIG. 3 gives by way of example the kinetics measured experimentally and the kinetics calculated by simulation for a few genes representative of the set of results obtained by the implementation of this variant.
  • Example 5 Modeling, simulations and validation of predictive capacity from a static graph of 133 molecules
  • This example shows the implementation of the whole method (steps A and B or A ', then C, D, E) and its predictive effectiveness in an application similar to the identification / selection of therapeutic targets.
  • a static graph corresponding to a network of molecular interactions in yeast was built according to the same principles as in Example 4. This graph more particularly includes the graph of Example 4, but with molecules and interactions additional. It includes 133 molecules, enzymes or transcription factors. It includes 407 uni- or bi-directional interactions between these molecules.
  • Glucose its concentrations over time were measured by the authors of the publication, at the same time as those for which the messenger RNA expression measurements were taken. The corresponding concentrations are given graphically in FIG. 4 of the cited article. In order to express these values as a ratio, each value of the concentration of Glucose in the culture medium at a given time was divided by the concentration of Glucose at the initial time of the experiment (this in order to measure the ratios by compared to the same frame of reference as for messenger RNA measurements, the ratios of which are expressed by the ratio of the measurement at time t divided by the measurement at the initial time).
  • Acetate and Glutamate the concentrations of these molecules were not measured by the authors. It was therefore decided to extrapolate values for these molecules from knowledge of the biological system studied and from the description of the experiments in the article. Insofar as this experiment is essentially based on the progressive fall in the concentration of glucose in the culture medium and since the other parameters of the culture medium are at first approximation considered to be constant, it was considered that the concentrations of Glutamate and of Acetate, respectively, were constant during the experiment.
  • Table 12 Values of the variables associated respectively with the molecule of the graph: Glutamate and with the molecule of the graph: Acetate
  • Steps A and B or step A ' were implemented in a similar manner to Example 4:
  • the initial state which corresponds to the experimental measurements at time TO, would be a stable state if the culture medium had been kept constant by supplementing it with glucose.
  • we did not define a final state of the graph which would correspond to a return to the initial state following a glucose supplementation after time 7 (TO + 12 hours).
  • the only experimental data to have been used for the calculation of the parameters were therefore the data actually described in the article and present on the website of Stanford University at the address: http://cmqm.stanford.edu /pbrown/explore/array.txt, and corresponding to the molecules of the network, without any extrapolation other than that concerning the molecules Glucose, Glutamate and Acetate described above.
  • the learning step for calculating the parameters was fixed at 16 minutes.
  • Step C was implemented as follows:
  • the messenger RNA expression screening data corresponding to this state were therefore those of the initial time used for the construction of the
  • Stage D was implemented as follows
  • This step consists in practicing iterative simulations, as described in the invention.
  • the yeast strain with the deletion of the TUP1 gene initially differs from the “wild” strain only in this deletion. This deletion returning to a constant inhibition of the expression of the TUP1 gene, the simulations consisted in simulating, in an iterative manner, the constant inhibition of each of the 133 molecules of the network, and in performing a calculation of propagation over the time of this inhibition within the network. For each simulation, a single molecule in the network was inhibited, since the state to be reached corresponds to an evolution of the modeled biological system (the yeast strain) following a single inhibition (deletion of the TUP1 gene). We therefore carried out 133 simulations.
  • the article's experimental data and the messenger RNA expression screening data concerning the strain exhibiting the TUP1 gene deletion (available on the Stanford University website at 1 'address: http://cmgm.stanford.edu/pbrown/explore/tupsearch.html).
  • the deletion of the gene was incomplete in this biological experiment, the ratio: [level of expression of the TUP1 gene in the deleted strain] / [level of expression of the TUP1 gene in the wild strain] being equal to 0.1 to a measurement , and 0.45 when replicating the measurement (average: 0.28). In the case of a complete deletion, this ratio would have been theoretically equal to 0, and equal in practice to the experimental measurement noise.
  • the inhibition was therefore defined numerically, for each molecule of the graph, as the multiplication by a factor of 0.1 of the level of expression of this molecule at the initial time (state to be modified as defined above), this factor corresponding to the value of the strongest inhibition measured experimentally for this gene.
  • the molecule i was different: the effect of each inhibition by a factor of 0.1 on each molecule in the graph was tested systematically.
  • the inhibition was imposed as a constant over time during the simulations: thus, during the propagation calculations, a possible propagation feedback on the inhibited molecule i
  • step D was implemented as described in the invention, without any notable particularity, and without presenting any particular difficulty for those skilled in the art.
  • the simulation calculations, consisting of propagating the initial inhibition over time, were performed using the same principles and the same tools as the simulations that are part of the parameter calculation procedure.
  • step D Each of the 133 simulations of step D thus resulted, at time 12 hours (duration of the simulated propagation), in the calculation of a new numerical value associated with each molecule of the network, defining a state of the graph: "state obtained by simulation ". We therefore obtained 133 different “states obtained by simulation”.
  • Step E was implemented as follows:
  • This step consists in prioritizing the molecules of the graph, and the effects exerted on these molecules during simulations, with reference to the greater or lesser proximity of the resultant of these effects with a state of the graph to be reached.
  • the state to be attained was the state of the yeast strain having a deletion of the TUP1 gene described above.
  • Table 14 Complete list of the values of the state to reach as defined above
  • Step E then consists in calculating the distance between on the one hand “the state to reach” of the graph, and on the other hand each of the 133 “states obtained by simulation” of the graph obtained in step D.
  • This distance calculation is described above (proximity of the states of the graph) and does not pose any particular difficulty for those skilled in the art. It consists in comparing two states of the graph by comparing two by two the set of values Xi associated with each molecule i of the graph.
  • Stage 1 consisted in carrying out a first classification by distance calculations in a conventional manner:
  • Order 1 distance sum of the absolute values of the differences between the values of Xj measured experimentally during the deletion of the TUP1 gene (X, 2 in the formula below) and the values of Xj measured by simulation (X ⁇ in the formula below) i
  • Stage 2 consisted, following this first classification, in refining it.
  • This classification step is given here by way of example but is not essential for the implementation of step E.
  • the sensitivity was then defined, for each of the 133 simulations, as the number of molecules of group Bj actually present in group A. This amounts to evaluating, for the quantitatively significant variations in expression of the molecules (greater than a factor of 2) if the simulation actually induces the variations present in the experimental data of the state to be reached. The higher the value of the sensitivity, the smaller the distance between the two states of the compared graph.
  • the specificity was then defined, for each of the 133 simulations, as the number of molecules of group Bi absent from group A. This amounts to assessing, for the quantitatively significant variations in expression of the molecules (greater than a factor of 2) if the simulation does not induce variations of expression absent in the experimental data of the state to be reached. The lower the value of the specificity, the smaller the distance between the two states of the compared graph.
  • the difference sensitivity - specificity therefore amounts to evaluating the distance on the combined criteria of induction by the simulation of the variations of expression present in the state to be reached and non-induction by the simulation of variations of expression absent in the state to reach.
  • sensitivity - specificity is also very simple and can for example be calculated manually, or automatic by an Excel software table (Microsoft) or any other equivalent software.
  • the relative overall error-learning calculation was 25.90%, which is satisfactory.
  • This global relative error result shows that the kinetics calculated during the simulations are close to the real data; random kinetics would have given an error calculation greater than 1.
  • FIG. 5 gives by way of example the kinetics measured experimentally (in white) and the kinetics calculated by simulation (in black) for a few representative molecules of the set of results obtained by the implementation of this variant of steps A and B or from step A '.
  • step B The calculation of the parameters carried out in step B therefore made it possible to obtain a dynamic graph which takes into account the experimental data used for the calculation.
  • the result of the first classification step is summarized by way of example in the following table, in the form of a summary table.
  • Each molecule in the network is classified by the distance between the target state of the graph and the state of the graph obtained by simulating constant inhibition by a factor 0.1 of this molecule.
  • Table 16 Distances between the target state of the graph and the state of the graph obtained by simulating constant inhibition by a factor 0.1 of each molecule
  • the classification in ascending order of the distances gives directly the classification of the molecules from the one whose inhibition is most likely to cause the state to reach of the graph to that whose inhibition is the less likely to cause it.
  • the TUP1 molecule is classified in first position, which is the expected result. The method is therefore validated.
  • FIG. 6 gives by way of example a schematic representation of this result of classification of the molecules of the network.
  • Each point corresponds to a molecule in the network.
  • the ordinates correspond to the calculated distance values.
  • On the abscissa the 133 molecules of the network are classified from left to right from that associated with the lowest distance to that associated with the highest distance.
  • the TUP1 molecule is classified in first position. This result is completely satisfactory.
  • the experimental test of the first 5 molecules as classified here would therefore give a success rate of 100% for the selection of the relevant molecule.
  • This classification also has the advantage of being able to define a “boundary” of all of the selected target molecules. Indeed, certain molecules of the graph do not send any interaction to another molecule (they are therefore “outputs” or “outputs” of the graph: “molecules outputs” in the rest of the text); the simulation of the inhibition of these molecules therefore does not cause the inhibition to propagate within the graph which therefore remains globally stable (since we considered that the initial state was stable). Therefore, these molecules are also classified in a contiguous group, from the 27th position to the 30th position. I_es molecules less well classified than the output molecules are therefore not of interest as a target molecule.
  • the classification can be interpreted as follows: the best-placed molecule is the one whose simulation of inhibition results in the state of the graph closest to the state to be reached. For the following molecules, we gradually move away from the state to reach, by going towards the state to modify that we reach when we arrive at the output molecules (this state not being modified during the simulation of the inhibition of output molecules, except for the output molecule). Beyond the output molecules, the inhibition simulations lead to states of the graph which gradually move away from both the state to be reached and the state to be modified. We did end up with the selection of a limited number of target molecules (those that are better classified than the outputs, here 26 molecules), which are themselves prioritized in terms of priority. The classification of the TUP1 molecule shows that this ranking is satisfactory.
  • Classification step 2 (the sensitivity - specificity calculation), although not essential given the above result, was then implemented in order to improve the classification of the targets. It was only applied to the molecules of the graph better classified than the output molecules in the previous step (26 molecules). The classification thus obtained is given in table 17.
  • Table 17 Final hierarchy of the selected target molecules (the 26 molecules before the outputs) The molecules are classified in descending order of the arithmetic difference [sensitivity - specificity]. We see that TUP1 is the only molecule in the graph for which the sensitivity is greater than the specificity. This individualizes it from others on qualitative criteria. Here again, TUP1 is in first position in the hierarchy of potential targets.
  • steps A and B, or A ', then C, D and E has been carried out more than ten times, and has systematically given similar results each time (with classification of TUP1 always in the first 5 molecules by the first classification mode). This shows without ambiguity the reproducibility of the method as well as its effectiveness.

Landscapes

  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Molecular Biology (AREA)
  • Physiology (AREA)
  • Biophysics (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Investigating Or Analysing Biological Materials (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Medicines That Contain Protein Lipid Enzymes And Other Medicines (AREA)

Abstract

La présente invention concerne le domaine de l'analyse intégrative des interactions moléculaires dans un système biologique. Elle porte en particulier sur un procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, permettant l'analyse desdites interactions lorsqu'un stimulus est appliqué au modèle dynamique, en vue notamment de hiérarchiser des molécules biologiques ou de sélectionner des cibles thérapeutiques vis-à-vis d'un problème biologique donné, pour en particulier définir une action thérapeutique à appliquer auxdites molécules. L'invention porte également sur un système informatique pour l'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, et l'analyse de ces interactions moléculaires lorsqu'un stimulus est appliqué au modèle dynamique, le système informatique comprenant au moins une unité centrale de traitement de données relié à au moins une base de données expérimentales quantitatives.

Description

METHODE ET SYSTEME DE SELECTION DE CIBLES THERAPEUTIQUES PAR L'UTILISATION DE RESEAUX DYNAMIQUES D'INTERACTIONS MOLÉCULAIRES
La présente invention concerne le domaine de l'analyse integrative des interactions moléculaires dans un système biologique. Elle porte en particulier sur des méthodes d'obtention et d'analyse des réseaux d'interactions moléculaires biologiques permettant, à partir de l'obtention de données expérimentales, d'identifier et de décrire les fonctionnements de ces interactions à la fois (i) entre des molécules interagissant deux à deux, (ii) au niveau des résultantes des interactions s'exerçant sur une molécule donnée, et (iii) au niveau de l'ensemble du réseau d'interactions considéré. Encore plus particulièrement, cette méthode d'analyse permet, une fois les fonctionnements de ces interactions décrits, de prédire les conséquences, sur l'ensemble du réseau d'interactions moléculaires considéré, d'actions d'activations ou d'inhibitions des molécules formant ce réseau. Elle permet ainsi notamment d'identifier des cibles thérapeutiques potentielles, de comprendre des mécanismes d'actions de xénobiotiques.
L'un des enjeux actuels majeurs des entreprises biotechnologiques et pharmaceutiques, est de développer de nouveaux médicaments, plus efficaces, contre la plupart des maladies, notamment (mais non uniquement) les maladies chroniques d'origines multi-factorielles qui représentent l'essentiel de la morbidité et de la mortalité dans les pays développés : maladies cardio-vasculaires, cancers, maladies psychiatriques et neurologiques, maladies immuno-allergiques, maladies endocriniennes (diabète...) etc. Les traitements actuellement disponibles pour la plupart de ces maladies ont des effets purement symptomatiques, souvent insuffisants même du point de vue purement symptomatique, sans action notable sur l'évolution propre de ces maladies, et avec souvent des effets indésirables importants. Par ailleurs, on ne dispose à ce jour d'aucun traitement réellement efficace pour certains syndromes ou maladies qui représentent des problèmes de santé majeurs, tels que les maladies neuro- dégénératives. La principale raison de cette situation est l'insuffisance actuelle de la compréhension des mécanismes physio-pathologiques aboutissant aux conditions pathologiques concernées, et notamment l'insuffisance de compréhension des mécanismes physiopathologiques moléculaires.
En effet, la plupart des médicaments existants ont été développés suivant une approche "pharmacologique" (désormais classique), consistant, schématiquement, à tester et sélectionner des molécules thérapeutiques potentielles (un grand nombre d'entre elles étant . obtenues par des méthodes de synthèse organique artificielle, notamment de type chimie combinatoire) sur des modèles physiopathologiques cellulaires et / ou animaux. Ces modèles sont censés reproduire tout ou partie des symptômes ou des modifications observées dans la pathologie. Cependant, une compréhension des mécanismes physiologiques, et notamment des mécanismes moléculaires, mis en jeu dans ces modèles, n'est pas requise pour leur mise en œuvre. Cette approche, fondée sur le crible à grande échelle de petites molécules de synthèse, présente donc l'avantage d'être relativement peu dépendante de la compréhension fine des processus physio-pathologiques impliqués dans les maladies concernées. Elle présente cependant une limite majeure, à laquelle elle est progressivement arrivée depuis environ deux décennies : sa dépendance vis à vis des modèles physiopathologiques utilisés, avec actuellement un épuisement des modèles génériques. Ceci est notamment lié au fait que la plupart de ces modèles ont été développés dans une logique d'interdépendance réciproque entre l'observation d'effets de molécules thérapeutiques et l'analyse progressive des actions pharmacologiques (moléculaires) de ces molécules. Ces modèles sont donc pour la plupart dépendants des effets pharmacologiques initialement observés, et ne permettent plus que de développer des médicaments d'effets proches de ceux déjà existants. Cette approche a progressivement évolué vers un coût élevé lié à un taux d'échec important dans le développement de nouveaux médicaments. Le développement de modèles animaux par transgénèse n'a pas, à ce jour, résolu ce problème de l'épuisement des modèles physiopathologiques : en effet, d'une part il s'avère que les gènes modifiés par transgénèse ne sont en général pas eux-mêmes des cibles thérapeutiques, et d'autre part l'approche de criblage de petites molécules de synthèse nécessite, pour être mise en œuvre de façon efficace, une orientation de la synthèse, soit par analogie avec des molécules existantes (ce qui ne permet le plus souvent pas d'innovation thérapeutique importante), soit par la connaissance préalable de la (ou des) molécule(s) cible(s), auxquelles on n'a pas directement accès par les modèles transgéniques. Par ailleurs, dans le cas des modèles d'animaux transgéniques de type knock-out, le fait qu'un gène cible thérapeutique éventuel ait été éliminé empêche tout criblage de molécules pharmacologiques potentiellement actives sur ce gène ou la protéine pour laquelle il code.
De ce fait, l'approche qui est de plus en plus mise en œuvre pour développer de nouveaux traitements pharmacologiques est une autre approche, dite "physiologique" ou "compréhensive", qui consiste à explorer et comprendre les mécanismes physiopathologiques, et notamment les mécanismes physiopathologiques moléculaires, aboutissant à la pathologie concernée, afin de définir les molécules de l'organisme à soigner, qui seront les molécules-cibles (ou "cibles thérapeutiques") des traitements chimiques. L'identification de ces molécules-cibles permet alors, dans un second temps, d'effectuer des cribles de molécules thérapeutiques potentielles de synthèse afin d'identifier celles qui vont modifier directement d'activité biologique de ces cibles thérapeutiques, ou encore d'effectuer des synthèses orientées de telles molécules thérapeutiques lorsque la structure spatiale des molécules-cibles est connue. Schématiquement, dans cette seconde approche, l'accent est mis sur la compréhension des mécanismes moléculaires physiologiques et physiopathologiques sous-tendant la maladie. Cette approche est elle aussi ancienne, puisqu'elle a débuté avec des techniques et méthodes de chimie organique permettant d'identifier et d'analyser des molécules intracellulaires (par exemple : description du cycle de Krebs), et elle s'est développée en intégrant les techniques et méthodes de la biologie moléculaire depuis environ 15 - 20 ans (séquençage d'acides nucléiques, clonage, transgénèse...). Il ne s'agit donc pas d'utiliser un modèle cellulaire ou animal physiopathologique pour cribler des petites molécules de synthèse à la recherche d'un effet thérapeutique, mais d'analyser d'abord le processus pathologique dans ces modèles (ou, lorsque c'est possible, directement chez l'Homme) afin de déterminer les cascades d'événements moléculaires menant à l'état pathologique pour identifier des cibles thérapeutiques potentielles. L'étape de synthèse de nouvelles molécules thérapeutiques n'intervient donc que dans un second temps. Cette approche "compréhensive" a aussi permis le développement de médicaments, mais elle s'est longtemps heurté à une limite majeure : jusqu'au milieu des années 1990, il était extrêmement difficile d'étudier plus d'une ou deux molécules en même temps. En conséquence, les cascades d'événements moléculaires décrites n'incluaient que quelques dizaines de molécules dans le meilleur des cas, alors que plusieurs dizaines de milliers de molécules différentes sont présentes dans une cellule eucaryote donnée. Cette approche ne permettait donc pas d'étudier les processus pathologiques en intégrant la complexité de ces processus, ceci particulièrement en ce qui concerne les maladies à déterminisme multi- factoriel. Enfin, sa logique tendait vers l'identification de cibles thérapeutiques uniques, mal adaptées pour soigner des maladies où une action sur une seule cible n'est pas suffisante. De fait, dans la plupart de ces maladies, la focalisation des traitements sur des actions pharmacologiques limitées s'est accompagnée de l'augmentation en miroir du nombre de médicaments prescrits à un même patient (ainsi, la plupart des maladies cardio-vasculaires, des maladies psychiatriques, etc. sont aujourd'hui traitées par poly-thérapie, et avec souvent des effets insuffisants et des interactions médicamenteuses difficiles à gérer et parfois néfastes). Les progrès enregistrés par cette approche ont essentiellement concerné la diminution des effets indésirables des traitements (ceux-ci ayant des actions moléculaires plus ciblées), sans être pour autant accompagnés d'une amélioration notable des effets thérapeutiques (un exemple en est le développement de la fluoxetine (Prozac ®). Dans ce contexte, il était crucial de pouvoir passer d'une analyse des processus pathologiques molécule par molécule à une analyse en parallèle de l'ensemble des molécules impliquées, permettant seule de rendre compte de la complexité de ces processus pathologiques. Ceci a été rendu partiellement possible depuis la fin des années 1990 par deux processus parallèles : d'une part l'identification d'une grande partie des molécules constitutives d'organismes vivants tels que l'homme et certains animaux modèles (séquençage de génomes entiers, identification en cours de l'ensemble des gènes présents dans ces génomes, déduction en cours des protéines correspondantes), et, d'autre part, le développement de techniques de biologie moléculaire qui permettent d'étudier un grand nombre de molécules différentes dans un même tissu. La plus significative de ces techniques, citée ici à titre d'exemple, est celle des puces à ADN permettant l'analyse en parallèle de virtuellement tous les ARN messagers présents dans un type cellulaire ou un tissu donné (Schena et al., 1995, Quantitative monitoring of gène expression patterns with a complementary DNA microarray, Science 270: 467-470; Lockhort et al., 1996, Expression monitoring by hybridization to high-density oligonucleotide arrays, Nature Biotechnology 14:1675-1680; Blanchard et al., 1996, Séquence to array: Probing the genome's secrets, Nature Biotechnology 14, 1649; brevet U.S. No. 5,569,588, publié le 29 Oct. 1996, Ashby et al., pour Methods for Drug Screening). Des méthodes d'analyse protéique à grande échelle ont aussi été développées, telles que l'utilisation d'hybrides dans la levure, le couplage électrophorèse 2D - spectrographie de masse, etc (McCormacket al., 1997, Direct analysis and identification of proteins in mixtures by LC,MS,MS and database searching at the low- femtomole level, Anal. Chem. 69(4):767-776; Chait, Trawling for proteins in the post-genome era, Nature Biotechnology 14:1544) ou sont actuellement en cours de développement (notamment la co-précipitation à grande échelle de protéines sur micro-colonnes).
Cependant, à ce jour, ces évolutions technologiques ont mené à la génération d'une grande masse de données biologiques, toujours croissante, sans que des techniques et méthodes satisfaisantes d'analyse et d'exploitation de ces données n'aient été développées. Ceci a entraîné le développement d'outils de biologie integrative pour les interpréter et sélectionner des cibles thérapeutiques pertinentes parmi une grande masse de données expérimentales générées.
Les méthodes de biologie integrative visent à analyser le rôle des molécules présentes dans l'organisme à soigner en tenant compte (donc en intégrant) dans cette analyse des autres molécules avec lesquelles elles interagissent. Leur objectif est donc de permettre d'obtenir des modèles des cascades (ou réseaux) d'interactions moléculaires du vivant, notamment de celles impliquées dans les processus pathologiques. Dans le contexte de la sélection de cibles thérapeutiques, de tels modèles visent à être appliqués pour sélectionner ces cibles. Plus précisément, une telle application doit permettre de prédire les conséquences des actions d'activation ou d'inhibition des molécules du réseau, afin d'identifier celles qui auront un effet thérapeutique. Il n'est envisageable de réaliser de telles prédictions à grande échelle et de façon suffisamment fiable que si le modèle permet de pratiquer des simulations systématiques des effets d'actions d'inhibition ou d'activation des molécules de la cascade. Les méthodes de modélisation proposées aujourd'hui sont d'une part, des méthodes produisant des modèles statiques et, d'autre part, des méthodes produisant des modèles dynamiques.
Les méthodes de modélisation produisant des modèles statiques consistent à construire des graphes statiques représentant des cascades d'interactions de molécules biologiques à partir de données de la littérature scientifique (publications dans des revues, analyse de profils d'expressions de molécules, prise en compte de données de séquences, etc.). Le graphe résultant peut être représenté sous la forme d'un schéma, le plus souvent en deux dimensions, dont les nœuds (ou sommets) du graphe sont les molécules, et où ces nœuds sont reliés par des trait ou des flèches (ou arcs, ou sommets du graphe) représentant les interactions entre les molécules. Des exemples de graphes statiques sont ceux construits dans diverses bases de données publiques telles que par exemple la base KEGG (M. Kanehisa and S. Goto : KEGG : Kyoto Encyclopedia of Gènes and Génomes. Nucleic Acids Research, 28(1) : 27-30, 2000). Cette méthode de modélisation aboutit à des résultats purement qualitatifs. Elle ne suffit pas à la mise en oeuvre de simulations quantitatives et dynamiques pour prédire les effets d'actions sur des cibles thérapeutiques potentielles. Cette limite est source d'un taux d'erreurs très important dans la sélection des cibles. De plus, il est extrêmement difficile pour un expert biologiste d'analyser de façon cohérente un graphe de plus de quelques dizaines de molécules, et cela devient impossible pour des graphes de plus d'une centaine de molécules. En conséquence, les cascades d'interactions moléculaires analysées sont de taille très réduite par rapport aux cascades réellement mises en jeu dans les organismes vivants, donc très incomplètes, et cette méthode ne permet pas de chercher des cibles de façon exhaustive. Mise en œuvre seule, elle est donc insuffisante au regard des enjeux cités plus haut. Dans les méthodes produisant des modèles dynamiques, les graphes statiques représentant les cascades d'interactions moléculaires sont utilisés pour créer des modèles dynamiques de ces graphes, reproduisant autant que faire se peut le comportement dynamique de la cascade biologique étudiée (ou voie biologique). Les méthodes utilisées à ce jour pour réaliser de tels modèles sont : Les méthodes qualitatives : - La méthode des réseaux booléens. - La méthode des formalismes logiques généralisés. - La méthode des formalismes fondés sur des règles (aussi appelés "rule-based' ou "knowledge-based'). Les méthodes probabilistes : - La méthode des équations stochastiques. - La méthode des réseaux Bayésiens. Les méthodes d'équations différentielles : - La méthode des équations différentielles ordinaires non linéaires. - La méthode des équations différentielles décomposées-linéaires (piecewise-linear differential équations). - La méthode des équations différentielles partielles et modèles de distribution spatiale.
Les méthodes mixtes : - La méthode des équations différentielles qualitatives.
Les principes sous-jacents à ces différentes méthodes sont résumés dans le tableau 1 ci-dessous.
Figure imgf000010_0001
Tableau 1 .Comparaison des méthodes de modélisation : principes sous-jacents
Ce tableau doit être lu en considérant les éléments suivants (1) Intégration de données quantitatives : certaines méthodes ne sont pas conçues pour utiliser et analyser des données quantitatives expérimentales biologiques (formalismes rule-based), ou les modifient de façon importante lorsqu'elles imposent une discrétisation des variables (réseaux booléens, etc.), d'où la notation : intégration "partielle". Ces méthodes ont toutes été initialement conçues pour s'affranchir au maximum de telles données. Ceci les limite dans leur fiabilité et dans leur possibilité d'application pour la recherche systématique de cibles thérapeutiques sur de grands réseaux.
(2) Formalisme : il s'agit des principes de représentation des interactions biologiques utilisés dans la méthode. On/off : les molécules sont soit présentes, soit absentes, sans état intermédiaire possible. Discrétisation des variables : le taux des molécules peut prendre un nombre limité de valeurs finies ; il s'agit d'un raffinement du formalisme précédent, mais qui représente mal la réalité biologique où les taux des molécules varient de façon continue. Probabilité de réaction chimique : spécifique des méthodes probabilistes où l'évolution du réseau est liée à la probabilité estimée des événements moléculaires individuels. Synthèse/ Dégradation : les effets des interactions sont représentées comme limités à des réactions de synthèse ou de dégradation des molécules, ces représentations étant celles des équations élémentaires de chimie, en général limitées à la loi d'action de masse (dont l'expression élémentaire est : si A+B- C, à l'équilibre : [C]=k1[A][B]). Diffusion : la diffusion des molécules dans le système biologique étudié ou hors du système biologique étudié (par exemple une cellule) est aussi prise en compte, comme équivalente à une synthèse ou à une dégradation (respectivement) au sein du système.
(3) Variables utilisées : toutes les méthodes existantes définissent les variables comme étant le taux, ou la concentration, ou la quantité totale, des molécules, noté ici Xj pour la molécule i , et non sa proportion de variation par rapport à un état étalon xi0. (4) Fonctions continues : pour une fonction continue, les variables changent de façon continue (comme c'est le cas dans les systèmes biologiques réels) et non discrète.
(5) Modèle déterministe : une fois le modèle calculé, le réseau ne peut passer d'un état à un autre que par un seul chemin (séquence unique d'états intermédiaires). Le fait qu'un modèle soit déterministe permet d'obtenir une croissance linéaire de la quantité de calculs lors des simulations. A l'inverse, dans les modèles non déterministes, la quantité de calcul requise lors des simulations tend à croître de façon exponentielle avec la taille du réseau, pouvant aboutir à une impossibilité de mise en œuvre pour de grands réseaux.
Du fait de leurs caractéristiques résumées dans le tableau 1 , ces différentes méthodes nécessitent des pré-requis et sont utilisables dans des applications qui sont résumées dans le tableau 2 ci-dessous :
Figure imgf000013_0001
antérieur : pré-requis et applications possibles
Le tableau 2 doit être lu en considérant les éléments suivants : (1) Connaissance fonctionnelle du réseau biologique :
Niveau A : connaissance de l'existence en soi des interactions moléculaires, et au moins une partie de leurs orientations et une partie des effets des interactions (activation/ inhibition ou synthèse/ dégradation). Seule la connaissance de niveau A est largement disponible à ce jour. Par conséquent, seule une méthode ne requérant qu'une connaissance de niveau A peut être appliquée à des réseaux étendus. Niveau B : niveau A avec toutes les orientations des interactions et de tous les effets des interactions. Niveau C : connaissance fonctionnelle étendue du réseau, c'est-à-dire niveau B plus d'autres données telles que : constantes des vitesses des réactions chimiques, description d'effets de seuil, description d'effets allostériques, etc. A ce jour, quel que soit l'organisme vivant considéré, pour la plupart des molécules des réseaux d'interactions moléculaires les connaissances de niveau C ne sont pas disponibles. Une description fonctionnelle détaillée du réseau biologique est nécessaire à la mise en œuvre de la méthode lorsque des connaissances de niveau C sont requises. Du fait de l'indisponibilité des connaissances de niveau C pour la plupart des molécules, toute méthode requérant ce type de connaissance pour sa mise en œuvre ne peut être appliquée qu'à de très petits réseaux bien étudiés et connus (quelques dizaines de molécules au maximum) et est de fait inadaptée à son application à de grands réseaux (de plus de 100 à 150 molécules).
(2) Puissance de calcul : Niveau A : croissance linéaire avec la taille du réseau (en nombre de molécules) de la quantité de calcul requise. Ceci correspond à la possibilité de mise en œuvre sur un serveur de puissance standard (grand public). Les méthodes mettant en œuvre des calculs dont la quantité croît de façon linéaire avec la taille du réseau peuvent être appliquées à des réseaux étendus (sous réserve de ne pas présenter d'autre limite à cette application). Niveau B : croissance de la quantité de calcul intermédiaire entre les cas A et C. Les méthodes mettant en œuvre des calculs dont la quantité croît de façon intermédiaire entre A et C sont théoriquement applicables à des réseaux étendus mais à un coût élevé voire très élevé (et sous réserve de ne pas présenter d'autre limite à cette application).
Niveau C : croissance exponentielle avec la taille du réseau (en nombre de molécules) de la quantité de calcul requise. Toute méthode mettant en œuvre des calculs dont la quantité croît de façon exponentielle avec la taille du réseau requiert une très grande puissance de calcul. A titre d'exemple, certaines applications des réseaux bayésiens nécessitent environ 30 minutes de temps de calcul sur un serveur équipé d'un processeur de 1 ,2 Giga Hertz pour un réseau de 8 molécules : sur un réseau de 32 molécules, le temps de calcul sur le même ordinateur serait dans ce cas de plus d'un an et demi. En pratique, même avec des ordinateurs les plus puissants actuels les méthodes présentant une croissance exponentielle du temps de calcul ne sont pas applicables à de grands ou très grands réseaux (quelques milliers à quelques dizaines de milliers de molécules et plus ; certaines d'entre elles ne sont pas applicables même à des réseaux de quelques centaines de molécules). (3) Taille maximale de réseau mis en œuvre : il s'agit de la taille maximale des réseaux sur lesquelles la méthode a pu être mise en œuvre à ce jour avec succès.
(4) Applicable à des réseaux de 1000 à 100000 molécules : cette possibilité d'application est liée (i) aux principes intrinsèques de la méthode (par exemple les réseaux Bayésiens, qui sont des réseaux linéaires et donc non adaptés à de grands réseaux biologiques comprenant des boucles de rétro- contrôle ne peuvent pas être appliqués à de grands réseaux), (ii) au niveau A, B ou C de connaissance fonctionnelle du réseau biologique requise, la nécessité d'une connaissance de niveau C rendant de fait la méthode inadaptée aux grands réseaux, et la nécessité d'une connaissance de niveau B la rendant très difficilement applicable à de tels réseaux, et (iii) à la puissance de calcul requise (niveaux A, B ou C), une croissance du temps de calcul de niveau C étant de fait non compatible avec une mise en œuvre sur de grands réseaux, et une croissance de niveau B rendant la méthode très difficilement applicable à de tels réseaux. (5) Mise en œuvre pour l'identification systématique de cibles thérapeutiques : il s'agit de la mise en œuvre effective de la méthode dans une recherche systématique de cibles au sein du réseau, sans a priori. Aucune des méthodes existantes n'a pu être mise en œuvre dans cette application à ce jour.
Toutes ces méthodes sont peu fiables dans leurs prédictions dès lors que le réseau dépasse une cinquantaine de molécules. Elles sont donc mal adaptées pour réaliser des modèles dynamiques corrects des réseaux d'interactions moléculaires des organismes vivants qui présentent les caractéristiques suivantes : - un grand nombre de types moléculaires différents sont impliqués : de quelques centaines à quelques dizaines de milliers, voire centaine de milliers, - les cascades mettent en jeu des boucles de rétro-action, avec une redondance des circuits, - les vitesses de propagation des activations / inhibitions des molécules au sein des réseaux sont différentes en fonction des circuits (i.e. des chemins de propagation au sein du réseau), - des réseaux extrêmement complexes difficiles à modéliser. Pour être réellement applicable à la prise en compte des données de génomique, transcriptomique et protéomique produites à grande échelle, dans un objectif d'identification systématique de cibles thérapeutiques, les modèles dynamiques construits doivent permettre de modéliser des cascades d'interactions moléculaires telles que décrites ci-dessus. Le fait de produire un modèle de la dynamique d'un réseau d'interactions moléculaires biologiques ne suffit pas en soi pour pouvoir sélectionner de façon fiable et rationnelle de nouvelles cibles thérapeutiques. A ce jour, toutes les méthodes développées n'ont pu être appliquées qu'à la simple description de processus moléculaires dans des réseaux biologiques de petite taille (quelques dizaines de molécules au plus) et à quelques simulations visant à reproduire des modifications connues du réseau. Aucune n'a été appliquée à la sélection systématique de cibles thérapeutiques parmi l'ensemble des molécules du réseau, y compris sur ces petits réseaux, et a fortiori sur de grands réseaux. En effet, une telle application requiert la mise en œuvre d'une stratégie de simulations appropriée, telle que décrite dans l'invention, et qui n'a pas été décrite avec les méthodes existantes (et pour certaines d'entre elles, n'est pas applicable même sur de petits réseaux). Cette application, à savoir la sélection de cibles thérapeutiques à partir de modélisations dynamiques des réseaux d'interactions moléculaires de grande taille effectivement mis en jeu dans les processus pathologiques n'est donc pas atteint par les méthodes décrites à ce jour. La présente invention a pour objectif de fournir une méthode d'obtention de modèles dynamiques de réseaux d'interactions moléculaires dans un système biologique, qui rendent possibles ce type d'applications. Pour une bonne intelligibilité de ce texte, un certain nombre de termes sont définis ci-dessous. Par interaction moléculaire entre deux (ou plus) molécules biologiques, il est entendu ici une interaction où une molécule (ou plus) active ou inhibe une autre molécule (ou plus). Le cas où une molécule d'un type donné interagit avec une autre molécule du même type n'est qu'un cas particulier de cette définition générale. Deux molécules sont définies ici comme étant du même type si elles ont la même formule chimique. L'activation (ou l'inhibition, respectivement) est définie ici comme l'augmentation (ou la diminution, respectivement) de l'activité biologique de la (ou des) molécule(s) sur laquelle (ou lesquelles) s'exerce l'interaction considérée. Cette augmentation (ou cette diminution, respectivement) de l'activité biologique peut correspondre soit à une augmentation (ou une diminution, respectivement) du nombre de molécules d'un type donné présentes dans le système biologique analysé, chacune gardant la même activité (ou fonction) biologique, soit à une augmentation (ou une diminution, respectivement) de l'activité des molécules d'un type donné, leur nombre restant constant, soit à une combinaison de ces deux mécanismes, soit à la résultante de ces deux mécanismes. L'activation (ou l'inhibition, respectivement) peut aussi être :1a conséquence d'une augmentation (ou d'une diminution, respectivement) du nombre de molécules associée à une diminution - (ou une augmentation, respectivement) de leur activité biologique, si la résultante globale en est une augmentation globale (ou une diminution globale, respectivement) de l'activité, et vice-versa.
L'activation (ou l'inhibition, respectivement) peut être non-nulle ou nulle en fonction des molécules considérées et du système biologique considéré. Elle peut être variable au cours du temps. Le fait que certaines interactions du réseau d'interactions moléculaires considéré correspondent à une activation (ou une inhibition, respectivement) nulle n'est qu'un cas particulier du champ de l'invention.
L'activité biologique d'une (ou de) molécule(s) biologique(s) considérée(s) correspond à toute capacité de la (ou des) molécules considérées à avoir une interaction chimique et/ou physique avec toute autre molécule d'un autre type (ou avec une autre molécule du même type). Cette interaction chimique et/ou physique peut résulter ou non dans l'acquisition (ou la perte) par une des molécules interagissant de capacités à avoir une interaction chimique et/ou physique avec toute autre molécule d'un autre type (ou avec une autre molécule du même type). Les interactions chimiques sont toute interaction entre deux molécules (ou plus) provoquant une réaction chimique (pouvant être représentée par une modification de la formule chimique d'une molécule, ou la synthèse, ou la dégradation d'une molécule). Les interactions physiques sont toute interaction entre deux molécules (ou plus) provoquant la formation d'un complexe stable ou instable entre ces molécules. Des exemples d'activités biologiques de molécules et d'interactions moléculaires correspondantes sont (de façon non exclusive) : l'activité d'activation de la transcription d'un gène donné (interaction moléculaire : protéine (facteur de transcription) - ADN), l'activité de mise en oeuvre d'une réaction chimique (interaction moléculaire : protéine (enzyme) - molécule (substrat), permettant la transformation de la molécule-substrat en molécule-produit de la réaction chimique), l'activité de formation d'un complexe moléculaire protéique ayant lui-même telle ou telle activité biologique (interaction moléculaire : protéine (sous-unité du complexe) - protéine (sous-unité du complexe)), etc. Par molécule biologique, il est entendu ici toute molécule, quelle que soit sa complexité, présente dans le système biologique considéré. Par système biologique, il est entendu ici tout organisme vivant, qu'il soit procaryote ou eucaryote, et qu'il soit unicellulaire ou pluricellulaire, et que le système biologique corresponde à cet organisme dans son entier ou à une partie de cet organisme. A titre d'exemples, on peut citer :
Organismes entiers - Une cellule (eucaryote ou procaryote) dans son ensemble. - Un ensemble de cellules interagissant directement ou indirectement entre elles, ou n'interagissant pas entre elles : 0 l'ensemble des cellules en culture dans une boite de Pétri ; 0 l'ensemble des cellules en formant un organe ou une partie de cet organe : noyau amygdalien d'un cerveau de mammifère. - Un être vivant pluricellulaire. - Les différents exemples plus leur environnement. Partie d'un organisme : - Un organelle d'une cellule, tel qu'une mitochondrie. - Un ensemble de molécules participant à une fonction biologique donnée, tel qu'un ensemble de molécules participant à la respiration cellulaire, ou un ensemble de molécules participant à la mort cellulaire, que cet ensemble de molécules soit constitué de toutes les molécules participant à ladite fonction biologique où une partie seulement d'entre elles.
L'ensemble des molécules formant le réseau d'interactions moléculaires tel qu'il est décrit sous forme d'un graphe statique dans la figure 2 est un exemple de système biologique. De nombreux graphes statiques sont par exemple disponibles dans la base de données publique KEGG (M. Kanehisa and S. Goto : KEGG : Kyoto Encyclopedia of Gènes and Génomes. Nucleic Acids Research, 28(1) : 27- 30, 2000). Tout système biologique est constitué de molécules, ces molécules interagissant les unes avec les autres de façon plus ou moins stable et variable au cours du temps et des effets de l'environnement de ce système sur le système biologique lui-même. A titre d'exemple, l'apoptose (mécanisme de mort cellulaire) est la résultante de l'interaction de multiples molécules (hormones, protéines, seconds messagers, etc..) qui, pour certaines d'entre elles, ont des interactions physiques ou chimiques plus ou moins stables au cours du temps.
Par réseau d'interactions moléculaires il est entendu ici l'ensemble des molécules analysées par la méthode de l'invention associé à l'ensemble (ou une partie de cet ensemble) de leurs interactions biologiques possibles. Le réseau peut comprendre toutes les molécules du système biologique concerné, ou seulement une partie de ces molécules. Pour plus de clarté, le réseau peut être représenté visuellement sous la forme d'un graphe (comme un exemple en est donné dans la description ci-dessous). C'est ce type de représentation visuelle qui est à l'origine de l'utilisation du terme de "réseau". Une telle représentation n'est cependant pas un pré-requis de l'invention. Le réseau peut aussi être représenté par un tableau (ou une matrice) dont par exemple chaque ligne correspond à une des molécules du réseau et dont les colonnes correspondent aux caractéristiques des interactions biologiques possibles de ces molécules (ou d'une partie de ces interactions ou de leurs caractéristiques). Un graphe est ici une représentation du réseau d'interactions moléculaires sous la forme d'un graphe dont les sommets (ou nœuds) correspondent aux molécules du réseau d'interactions moléculaires représenté et dont les arrêtes (ou arcs) reliant les sommets correspondent aux interactions moléculaires du réseau d'interactions moléculaires représenté. Dans la suite du texte, il sera très souvent fait référence à un tel graphe, bien qu'il ne soit pas indispensable d'en réaliser un physiquement. Etant donné qu'il ne s'agit que d'une représentation symbolique du réseau, une référence au graphe correspond en réalité à une référence au réseau. Par variable associée à un sommet du graphe, il est entendu ici une variable quantitative au sens mathématique du terme, pouvant prendre des valeurs numériques, et dont la valeur à un état donné du graphe représente l'état du sommet correspondant en ce qui concerne une quantité se rapportant à une molécule du système biologique considéré. Suivant les cas, cette quantité peut être un niveau d'expression d'un gène exprimé dans le système biologique (par exemple, l'abondance d'ARN messagers, mesurable notamment par la technique des puces à ADN), un niveau d'abondance d'une protéine, un niveau d'activité d'une protéine, un niveau d'abondance d'un métabolite, etc, pourvu que la quantité considérée soit mesurable expérimentalement, par un moyen direct ou non. Un état d'un graphe est un graphe pour lequel une valeur numérique est donnée pour chaque variable (associée à chaque sommet). Le cas où une valeur numérique non nulle n'est donnée que pour une partie des variables (et associée aux sommets correspondants), une autre partie des variables (associées à d'autres sommets) étant nulles, n'est qu'un cas particulier d'état du graphe. Un état du graphe donné est une représentation d'un état réel ou simulé du réseau d'interactions moléculaires correspondant, et par extension une représentation d'un état réel ou simulé du système biologique correspondant. A titre d'exemple, dans certaines représentations d'un réseau d'interactions moléculaires sous la forme d'un graphe, le fait de donner à une variable associée à un sommet du graphe la valeur nulle peut correspondre à une représentation de la situation où la molécule correspondant à ce sommet n'est pas présente dans le réseau d'interactions (ce qui ne signifie pas qu'elle n'est pas présente dans le système biologique), ou bien de la situation où son activité biologique est nulle. Le fait de donner une valeur nulle à un certain nombre de variables correspond donc à considérer qu'à un temps donné celles-ci n'interagissent pas avec le reste du réseau, mais leur valeur peut devenir non nulle à un autre temps suite à une modification de l'état du réseau. Le fait de donner une valeur nulle à une variable ne revient donc pas nécessairement à exclure le sommet correspondant du réseau. Dans certains cas particuliers, il est possible de donner une valeur constamment nulle à un certain nombre de variables, ce qui correspond alors à exclure les sommets correspondants du réseau et donc à travailler sur un sous-réseau. Pour travailler sur un sous-réseau, on préférera toutefois faire une hypothèse conservatrice, c'est-à-dire considérer la valeur des variables exclues comme constante, ce qui permet de ne pas modifier la structure du réseau.
L'invention concerne également un système informatique pour l'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, et l'analyse de ces interactions moléculaires lorsqu'un stimulus est appliqué au modèle dynamique, comprenant au moins une unité centrale de traitement de données reliée à au moins une base de données expérimentales quantitatives, le système informatique comprenant : A) un module de construction d'un graphe statique, dont les sommets représentent des molécules biologiques et les arcs représentent des interactions physico-chimiques existant entre ces molécules, chaque sommet étant associé à une variable quantitative mesurée expérimentalement et chaque arc du graphe étant associé à une relation mathématique; et B) un module d'apprentissage pour calculer les paramètres de chaque relation à partir des données expérimentales quantitatives concernant les sommets du graphe, par la mise en œuvre de techniques d'apprentissage par descente de gradient utilisées pour le paramétrage de réseaux.
Le système informatique selon l'invention peut en outre comprendre : C) un module de simulation pour effectuer plusieurs procédures itératives de simulation consistant à imposer un stimulus à un état de graphe mesuré expérimentalement et choisi comme « état à modifier », le stimulus modifiant la valeur d'une ou de plusieurs des variables quantitatives associées aux sommets du graphe, constituant ainsi un état de départ de la simulation à partir duquel un calcul de propagation est effectué au sein du graphe, pour l'obtention d'un « état final du graphe »; et D) un module d'itération pour la modification du stimulus.
Le système informatique selon l'invention peut en outre comprendre : E) un module de calcul de proximité entre I' « état final d'un graphe » et I' « état à modifier », ou entre I' « état final d'un graphe » et un état voulu, et de hiérarchisation des sommets et des stimuli imposés sur les sommets du graphe, les sommets hiérarchisés correspondant à des cibles thérapeutiques classées.
Le système informatique selon l'invention forme un outil d'analyse de données expérimentales biologiques, et notamment un outil de hiérarchisation de molécules biologiques vis-à-vis d'un problème biologique.
La présente invention a entre autres pour objet d'apporter des solutions techniques aux difficultés exposées plus haut, notamment en apportant la possibilité de construire des modèles dynamiques utilisables pour des réseaux d'interactions moléculaires de plus de 100, de plus de 200 molécules ou même davantage, dans les applications décrites. Un premier aspect de l'invention est un procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires . dans un système biologique, permettant l'analyse desdites interactions et plus précisément permettant l'analyse dudit réseau d'interactions lorsqu'un stimulus est appliqué au modèle dynamique, en vue notamment de hiérarchiser des molécules biologiques ou de sélectionner des cibles thérapeutiques vis-à- vis d'un problème biologique donné, pour en particulier définir une action thérapeutique à appliquer auxdites molécules, ledit procédé étant mis en œuvre par un système informatique et comprenant les étapes suivantes : A) à partir d'un graphe statique dont les sommets représentent des molécules biologiques et les arcs représentent des interactions physico-chimiques existant entre ces molécules, associer une variable quantitative Xi mesurée expérimentalement à chaque sommet i, et une relation mathématique à chaque arc du graphe, chacune desdites relations présentant les caractéristiques suivantes : - elle comprend un terme inertiel (i) qui tend vers une limite finie; - elle comprend un terme (ii) tendant à faire revenir les variables Xj à leur état initial, de signe inverse au terme inertiel (i), et dont la variation en fonction du temps croit en valeur absolue de façon plus lente que la variation en fonction du temps du terme inertiel (i); - elle comporte un facteur de pondération wy qui permet de tenir compte de la combinaison d'effets pouvant s'exercer sur chaque sommet du graphe; B) calculer les paramètres de chaque relation à partir de données expérimentales quantitatives concernant les sommets du graphe, par la mise en œuvre de techniques d'apprentissage par descente de gradient utilisées pour le paramétrage de réseaux.
Le signe réel du terme (ii) est déterminé par le résultat du calcul de son ou ses paramètre(s). Ce terme (ii) est de signe inverse au terme (i) une fois les paramètres calculés, mais cela n'apparaît pas obligatoirement dans sa formulation mathématique, où l'on ne précise pas a priori le signe du ou des paramètre(s) associé(s).
Dans une mise en œuvre préférée du procédé ci-dessus, chaque variable quantitative associée à un sommet représente la variation relative de la quantité de la molécule correspondant audit sommet, par rapport à un état étalon du système biologique. Comme mentionné ci-dessus, la "quantité de la molécule associée à un sommet" peut concerner n'importe quel aspect mesurable directement ou non ce cette molécule, qu'il s'agisse de sa concentration, son activité, son taux d'expression, etc. Dans cette variante où les Xj sont des rapports à un état étalon, ledit état étalon est de préférence un état stable du système biologique, dans lequel la quantité de chaque molécule associée à un sommet du graphe est mesurable expérimentalement. Comme reprécisé dans la description d'une mise en œuvre pratique ci-dessous, cet état étalon peut correspondre à un état physiologique donné (par exemple sain ou malade) réellement observable, ou à un état artificiel du système, par exemple à l'état d'un pool de plusieurs échantillons biologiques prélevés dans des conditions expérimentales différentes. Les variations relatives de quantité des molécules du réseau sont donc représentées sous la forme de variables dépendantes des variations relatives de quantité des molécules interagissant sur elles (i.e. en interaction avec elles et en amont dans le réseau en termes de propagation des activations / inhibitions). La définition des variables correspond directement aux mesures expérimentales disponibles : en effet, dans la plupart des technologies de biologie moléculaire (dont les criblages d'expression d'ARN messagers), la quantité absolue des molécules présentes dans le système biologique d'intérêt .-n'est pas mesurée (ni mesurable) ; seule la proportion de leur variation par rapport à un état de référence est mesurable.
Soit les n molécules j (1-»n), représentée par les n sommets j (1-»n) du réseau, interagissant sur la molécule i, représentée par le sommet i du réseau. Dans les procédés d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon l'invention, les termes inertiel (i) et de retour à l'état initial (ii) permettent de calculer les valeurs de Xi et les variations de valeurs de Xi au cours du temps en fonction des valeurs des Xj (1-»n) et des variations des valeurs des Xj (1-»n) au cours du temps.
Par l'expression "terme inertiel", on entend : - une résistance au changement, résistance notamment initiale, et - un délai pour arriver à la variation maximale, ce qui permet de rendre compte des complexités des propagations dans le réseau.
En particulier, le terme inertiel (i) a pour objet de permettre d'intégrer une résistance des variables au changement et un décalage temporel entre les modifications des variables en amont et aval du réseau. Il introduit en particulier : - l'intégration du facteur temps - la prise en compte des différences de vitesses de propagation au sein du réseau en fonction des sous-circuits - la prise en compte des retards temporels consécutifs aux influences des boucles de rétro-contrôle sur la propagation dans le réseau, et - il permet de calculer les cinétiques des interactions moléculaires au sein du réseau directement à partir des données expérimentales, sans connaissance préalable des constantes de vitesse de ces cinétiques, et sans faire d'à priori sur d'éventuels autres paramètres.
Ce terme inertiel (i) tend vers une limite finie, ce qui permet d'éviter des divergences importantes lors des simulations (amélioration de leur fiabilité) : ceci évite le risque de divergence (ou d' "explosion") des valeurs des variables liées à des propagations itératives dans des boucles de rétroaction ou lors de simulations portant sur des temps prolongés. Le fait de pouvoir, en évitant de telles divergences, obtenir des convergences satisfaisantes lors de simulations portant sur des durées longues (qui soient en rapport par exemple avec les temps d'installation de processus pathologiques), est une caractéristique importante de l'invention. La formulation de ce terme inertiel est de préférence peu contraignante, et permet de rendre compte de formes multiples de relations. Pour cela, il peut être avantageusement exprimé sous la forme d'une relation mathématique présentant une ou plusieurs inflexion(s), ce qui permet de limiter les contraintes imposées aux modèles et de pouvoir pratiquer des modélisations fiables dans les situations où la forme des cinétiques n'est pas connue a priori, ce qui est une situation constante dès que l'on modélise un grand réseau (plus d'une centaine de molécules). Des exemples de telles sous-relations mathématiques pouvant être utilisées sont les relations sigmoïdes, les relations d'oscillation, et, d'une façon générale, toute fonction mathématique tendant à une ou des limite(s) finie(s) et pouvant être infléchie. Le terme (ii) tendant à faire revenir les variables à leur état initial (ou d'équilibre antérieur), permet de rendre compte des phénomènes d'homéostasie et de l'existence d'états d'équilibre du réseau, tout en diminuant de façon significative les risques de divergence lors des simulations (amélioration de leur fiabilité). Une fois les paramètres des relations mathématiques calculés, il est de signe réel inverse au terme inertiel (i), et sa variation au cours du temps croit en valeur absolue de façon plus lente (Le., de façon temporellement plus tardive) que la variation en fonction du temps du terme inertiel (i). Par le terme (i), Xj et les variations de Xj dépendent des Xj (1→n) et des variations des Xj (1-»n). Le terme (i), qui fait tendre Xj vers une valeur finie, est donc exprimé en fonction des Xj (1-»n).
Le terme (ii) est, lui, exprimé en fonction de Xj (et non des Xj (1-»n)). La valeur de ce terme ne peut donc changer que si la valeur de Xj change, celle-ci changeant si les valeurs des Xj (1-»n) changent.
Toute variation initiale de l'effet du terme (ii) sur la valeur calculée de Xj peut donc être considérée comme consécutive à une variation préalable de l'effet du terme (i) sur la valeur calculée de Xj. Ceci s'applique notamment si l'on considère qu'il existe un état stable du réseau ; à l'état stable, les termes (i) et (ii) s'équilibrent, de telle sorte que Xj reste constant ; à partir de cet état, toute variation de Xi est consécutive à une situation où l'effet du terme (i) sur la variation de Xj est plus grand en valeur absolue que l'effet du terme (ii) sur la variation de Xj. En effet, une fois les paramètres des termes (i) et (ii) calculés, le terme (ii) calculé est de signe opposé au terme (i) calculé, et tend, lors du calcul des valeurs de Xj à diminuer l'effet du terme (i) sur les variations des valeurs de
Xi.
Par conséquent, Xj ne peut présenter une variation que si, à un temps donné au moins, la variation de Xj au temps suivant calculée par le terme (ii) est inférieure en valeur absolue à la variation de Xj au temps suivant calculée par le terme (i). En d'autres termes, Xj ne peut présenter une variation, à partir d'un état stable, que si, sur un espace de temps donné au moins, la variation de la valeur calculée du terme (ii) est inférieure en valeur absolue à la variation de la valeur calculée du terme (i). Cette caractéristique est inhérente au fait que le terme (i) est exprimé en fonction des Xj (1→n) alors que le terme (ii) est exprimé en fonction de Xj. A partir d'un état stable, la variation du terme (ii) est initialement inférieure en valeur absolue à la variation du terme (i). Lors de l'évolution de la valeur de Xj au cours du temps, l'effet du terme (ii) sur la variation de Xj peut, ou non, devenir supérieur à l'effet du terme (i) sur la variation de Xj. Si c'est le cas, Xj va tendre à retourner vers sa valeur initiale.
En fonction des valeurs des paramètres calculés des termes (i) et (ii), des valeurs des Xj (1→n) et des valeurs de Xj, Xj peut éventuellement retourner à sa valeur initiale, ceci notamment si les Xj (1→n) retournent à leur valeur initiale.
Si un stimulus est appliqué de façon constante sur un ou plusieurs sommets du réseau, on peut cependant aboutir à une situation où les Xj (1-»n) ne retournent pas à leur valeur initiale. Dans ce cas, Xj peut ne pas retourner à sa valeur initiale. Si, à un temps donné, les effets des termes (i) et (ii) sur la variation de Xj s'équilibrent à nouveau, on aboutira alors à une nouvelle stabilité de Xj, à une valeur différente de sa valeur initiale. La méthode permet donc de rendre compte du passage du réseau d'un état stable donné à un autre état stable, différent. Elle permet aussi de rendre compte de l'évolution du réseau lors d'états instables.
Enfin, comme le terme (i) fait tendre Xj vers une limite finie, et comme le terme (ii) est exprimé en fonction de Xj, le terme (ii) est contraint par Xj : par la résultante des termes (i) et (ii) la valeur calculée de Xj ne peut sortir d'un intervalle fini. Cette caractéristique (Xj tendant vers une limite finie par le terme (i) et expression du terme (ii) en fonction de Xj) permet de rendre compte d'états stables, et de contraindre les valeurs de Xj dans un intervalle fini.
Le fait de pouvoir calculer les paramètres des relations associées aux arcs du graphe directement à partir des données expérimentales, sans nécessiter d'hypothèse préalable ou de fixation à des valeurs arbitraires, est rendu possible par l'utilisation de relations de forme peu contraignante, ne requérant pas une connaissance préalable des cinétiques des interactions moléculaires.
Comme mentionné ci-dessus, les procédés d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon l'invention comportent une deuxième étape (étape B), dans laquelle on calcule les paramètres des relations associées à chacun des arcs du graphe, à partir de données expérimentales quantitatives concernant les sommets du graphe. Ce calcul est effectué de préférence par la mise en œuvre de techniques d'apprentissage. On obtient alors un graphe dynamique, entièrement déterministe, consistant au graphe statique aux arrêtes duquel sont désormais associées des relations mathématiques dont les paramètres ont tous été définis numériquement.
Cette étape de calcul peut être effectuée par l'utilisation de procédures d'apprentissage utilisées pour le paramétrage de réseaux en intelligence artificielle, par exemple celles développées en informatique dans les méthodes de "réseaux de neurones" (dont les réseaux de neurones récurrents) par descente de gradient "simple" (en prenant comme base de calcul les couples de données (Xj, Xj) fournis par les données expérimentales indépendamment les uns des autres), ou par descente de gradient dans le temps (où ces couples ne sont pas considérés comme indépendants). Les couples de données (Xj, Xj) fournis par les données expérimentales sont définis comme suit : soit i une molécule du réseau, représentée par le commet i, et soit j toute molécule du réseau interagissant sur i, représentée par le sommet j. Xj et Xj sont les variables associées aux sommets i et j, respectivement. Les mesures expérimentales des valeurs des Xj et des Xj dans des conditions expérimentales données et à des temps expérimentaux donnés permettent d'obtenir des valeurs numériques des Xj et des Xj. Un couple de données expérimentales (Xj, Xj) correspond aux valeurs mesurées de Xj et Xj à un état expérimental donné (même temps, même condition expérimentale).
Les données expérimentales utilisées pour réaliser l'étape B) mentionnée ci-dessus présentent les caractéristiques suivantes : Nature des données expérimentales. Ces données . sont des données quantitatives concernant les molécules (correspondant aux sommets du graphe) et sont par exemple des niveaux d'expression de gènes exprimés dans le système biologique (par la mesure de l'abondance d'ARN messagers, par exemple par la technique des puces à ADN) et / ou des niveaux d'abondance de protéines et /ou des niveaux d'activité des protéines et / ou des niveaux d'abondance de métabolites. Comme précisé plus haut, ces données sont exprimables sous la forme d'une proportion de variation de quantité par rapport à une situation de référence (état étalon).
Compilation des données de réseau statique (ou graphe statique : identification d'interactions j → i) et des données expérimentales (mesures de valeurs des variables Xj). Ces données peuvent être extraites de la littérature scientifique au sens large, ceci incluant les bases de données biologiques publiques ou privées (telles que par exemple la base de données "TRANSFAC du "German Research Centre for Biotechnology "
(GBF) accessible par l'adresse internet : http://transfac.qbf.de/ (Wingender et al., 2001 , The TRANSFAC System on gène expression régulation, Nucleic Acids Research, 29 (1) : 281-283), ou encore la base de données " BIOMOLECULAR INTERACTION NETWORK DATABASE " (BIND) de l'université de Toronto, accessible par l'adresse internet : http://www.bind.ca (Bader et al., 2003, BIND : the biomolecular interaction network database, Nucleic Acids Research 31 : 248-250), ou encore la base de données KEGG (M. Kanehisa and S. Goto : KEGG : Kyoto Encyclopedia of Gènes and Génomes. Nucleic Acids Research, 28(1) : 27- 30, 2000), ou bien être générées par des expériences de biologie moléculaire dédiées, notamment par l'utilisation des techniques de criblages à grande échelle. En fonction du système biologique d'intérêt, des molécules formant le réseau d'interactions moléculaires, de la problématique scientifique biologique (étude d'un modèle de maladie, étude de toxicité d'un produit, étude de processus du développement, etc.), des paradigmes expérimentaux adéquats ou disponibles (cultures de cellules, étude de tissus, etc.), la personne de l'art définira le type de données expérimentales d'intérêt. Des exemples du type de données utilisables en fonction des applications de l'invention sont donnés ci-après, dans la description de la méthode des simulations. L'homme du métier mettra donc en œuvre la ou les méthodes de compilation lui convenant le mieux pour effectuer cette étape, qui intervient en amont de la méthode d'analyse constituant la présente invention.
Enregistrement des données expérimentales. Ces données expérimentales sont enregistrées avantageusement dans une base de données, de nombreux systèmes de bases de données existant pour ce faire et pouvant être mis en œuvre et utilisés de façon simple par toute personne de l'art du domaine de la bio-informatique (bases de données commerciales Oracle, Microsoft SQL server, FileMaker, bases de données d'accès libre : postgreSQL). Ces données peuvent aussi être enregistrées sous le format d'un tableau ou d'un fichier plat.
Indexation des données expérimentales. Ces données expérimentales peuvent être indexées automatiquement au graphe. Le rôle de cette indexation est de relier chaque donnée expérimentale à l'objet biologique correspondant du graphe (sommet du graphe, ou arrête du graphe pour les couples de données (Xj, Xj), de façon à pouvoir utiliser conjointement ces deux types d'informations (données expérimentales et graphe) lors de la mise en œuvre du système de calcul des paramètres. De nombreux systèmes de bases de données commerciaux ou gratuits permettent de créer cet indexage sans difficulté technique particulière pour l'homme de l'art du domaine de la biologie ou de la bio-informatique (bases de données commerciales Oracle, Microsoft SQL server, FileMaker, bases de données d'accès libre : postgreSQL). Alternativement, si les données concernant le graphe et les résultats expérimentaux ont été enregistrées sous le format de tableaux ou de fichiers plats, ou d'un tableau ou fichier plat commun, ces données étant liées de fait dans ce cas, cette étape d'indexation peut ne pas être nécessaire en soi.
Forme des données expérimentales. Dans une mise en œuvre préférentielle, les données expérimentales des valeurs des couples (Xj, Xj) sont sous la forme de cinétiques d'expression. Par cinétique d'expression il est entendu ici un ensemble de séries de données expérimentales ordonnées dans le temps, chaque série de données correspondant à un ensemble de valeurs de couples (Xj, Xj) mesurés expérimentalement à un temps donné. Chaque série de données peut concerner soit l'ensemble des sommets du graphe, soit uniquement un sous-ensemble de ces sommets. Les différents temps correspondent à des temps successifs au cours de l'observation d'un processus biologique mettant en oeuvre le système biologique modélisé par le graphe, que ce processus soit naturel ou induit artificiellement en laboratoire. Une telle cinétique comprend de préférence au moins trois temps successifs, et, pour améliorer la qualité du calcul des paramètres, plus de trois temps.
Plusieurs cinétiques indépendantes, correspondant à des processus biologiques différents (Le., mettant en jeu des sous-réseaux différents d'un même réseau global, ces sous-réseaux pouvant ou non présenter des parties communes), peuvent être utilisées simultanément. Ceci peut permettre d'améliorer la qualité du calcul des paramètres, et donc la qualité des simulations.
Dès lors qu'au moins une cinétique d'expression est disponible, il est possible d'utiliser simultanément aussi des données expérimentales des valeurs des couples (Xj, Xj) obtenues par des expériences indépendantes les unes des autres (sans description de cinétiques d'évolution du système biologique étudié au cours du temps).
La méthode de calcul des paramètres des relations, à l'étape B) des méthodes de l'invention, tient de préférence compte des principes suivants :
Mesure expérimentale d'un état stable du système biologique. Le graphe est considéré comme étant dans un état stable de référence à un temps donné, cet état stable étant mesurable expérimentalement. L'état stable de référence en question correspond à un état existant et mesurable du système biologique étudié, pouvant être considéré comme stable dans le temps vis-à-vis du processus biologique modélisé. Bien qu'un système biologique soit le plus souvent, du fait de ses interactions avec l'environnement et de ses rythmes biologiques propres, en train de se modifier, on peut définir, du fait de l'existence des processus homéostatiques, des états où ces modifications sont au maximum "oscillantes" autour d'états homéostatiques, et a priori de faible amplitude. Dans cet état, le processus modélisé n'est pas lui-même en train d'évoluer significativement.
Cet état ne doit pas être confondu avec l'état étalon. L'état étalon, qui est défini arbitrairement par l'expérimentateur biologiste sert à effectuer des mesures quantitatives expérimentales. L'état stable de référence correspond à un état réel du système modélisé (Le., non artificiel), et sert de référence pour le calcul des paramètres du modèle. Il est considéré comme un état du système où les processus d'activation et d'inhibition au sein du réseau sont équilibrés, ou présentent des oscillations faibles autour d'un état d'équilibre théorique. Il représente l'état vers lequel le système tend en général à retourner lors des simulations. Il peut être le même, ou différent, de l'état étalon. L'état stable de référence est directement mesurable expérimentalement dès lors que le problème biologique étudié permet de définir un état de référence du système biologique.
A titre d'exemple, une culture cellulaire dont le nombre de cellules est arrivé à un plateau (absence de divisions cellulaires) et dans un milieu de culture stable, avant toute induction de stimulus, ou un animal adulte sain avant toute induction de processus pathologique, peuvent être considérés comme des états stables de référence. Dans le premier cas, les cascades d'interactions moléculaires mises en jeu par le stimulus dont on cherche à modéliser les conséquences ne sont pas activées au-delà des processus homéostatiques. Dans le second cas, les cascades mises en jeu par le processus pathologique à modéliser ne sont pas non plus en œuvre : l'état de référence est stable vis à vis du processus biologique modélisé. L'état stable ne doit pas nécessairement être l'état initial du système biologique dans le cadre du processus biologique étudié. Dans un autre exemple, l'état sain peut être considéré comme un état stable initial de référence si l'on étudie l'installation d'un processus pathologique à partir de cet état sain.
La mesure des Xj de l'ensemble des sommets du graphe dans cet état est utilisée, dans le calcul des paramètres, comme référence stable du graphe, notamment pour la procédure de minimisation des erreurs.
L'état stable est défini mathématiquement par le vecteur de l'ensemble des valeurs expérimentales des variables de chaque sommet mesurées à l'état biologique correspondant (mesures effectuées pour tous les sommets du graphe). Dans une mise en œuvre préférentielle, l'état étalon pour les mesures est l'état stable. Dans ce cas, comme les variables sont définies par (voir l'exemple 1 de mise en œuvre) : Xj = Xjt/Xjo, en théorie, à l'état stable, puisque le réseau ne se modifie pas, quel que soit t, x,t = XJO, donc Xj = 1 , pour tout sommet i. C'est le fait d'induire une modification du réseau par l'application de stimuli lors des expériences biologiques qui va "déstabiliser" le réseau, aboutissant à la mesure de cinétiques où Xu ≠ XJO et Xj ≠ 1.
Dans cette mise en œuvre, on peut donc éventuellement définir un état stable arbitraire où quel que soit i, Xj = 1.
En pratique, lors de la mesure expérimentale de cinétiques, au premier temps (t0), les Xj sont proches de 1 en général (si l'état étalon de mesure est le temps t0).
Mais ce qui est important n'est pas tant le fait que les Xj soient égaux à 1 en théorie et proches de 1 lors des mesures expérimentales, mais le fait en soit que cet état soit considéré comme stable. En effet, lors du calcul de paramètres des relations mathématiques entre les Xj et les Xj par des techniques de réseaux de neurones avec minimisation des erreurs, le fait de définir un état comme stable (au moins au début de la cinétique) introduit une contrainte forte dans le calcul des paramètres et améliore ainsi significativement leur calcul. Pour que le modèle obtenu soit pertinent vis-à-vis du ou des processus biologique(s) étudié(s), il est préférable de s'assurer que cet état stable existe biologiquement, en le validant, par sa mesure expérimentale. Si l'état stable est différent de l'état étalon, les valeurs des Xj à l'état stable ne peuvent être définies rationnellement que par leur mesure expérimentale. Il est également possible de décider arbitrairement de le définir par V i, X, = 1 , et d'introduire (au sens "ajouter") ce vecteur des Xj au temps initial des cinétiques sans l'avoir mesuré. Ceci revient à considérer l'état étalon comme stable, arbitrairement. Ceci est souvent possible si l'état étalon ne correspond pas à un pool de tissus biologiques différents. Dans une mise en oeuvre préférentielle, les données expérimentales sont mesurées au cours d'une cinétique (voir plus haut). Dans le cas où le processus biologique d'intérêt est étudié au cours du passage d'un état stable initial à un état stable final, et où des mesures expérimentales sont effectuées à ces deux états et à des temps intermédiaires, deux états stables sont définis : l'état initial et l'état final de la cinétique du processus biologique étudié. Cependant, le fait de disposer de mesures expérimentales correspondant à deux états stables n'est pas un pré-requis à la mise en oeuvre de l'invention.
Le fait de définir un état stable n'est pas non plus un pré-requis à la mise en œuvre de l'invention.
Lissage des données : si l'ensemble des données expérimentales est très restreint, une procédure de lissage des données expérimentales peut être mise en œuvre préalablement au calcul des paramètres, pour permettre d'augmenter le nombre de valeurs de couples (Xi; Xj) disponibles, en calculant des valeurs intermédiaires de ces couples à partir de la courbe lissée. Cette procédure, classique, ne pose pas de difficulté particulière à l'homme de l'art.
Le calcul des paramètres des relations (Xj, Xj) est effectué par la mise en œuvre de techniques d'apprentissage utilisées pour le paramétrage de réseaux en intelligence artificielle (telles que celles mises en oeuvre pour les réseaux de neurones), à partir des données expérimentales quantitatives concernant les variables du graphe.
A titre d'exemple, ce calcul peut utiliser des algorithmes de résolution numérique de propagation ou de rétro-propagation avec calcul de l'erreur.
Des paramètres sont arbitrairement fixés, une propagation ou une rétro- propagation est effectuée, puis l'erreur est calculée entre les résultats calculés et les résultats expérimentaux. Les paramètres sont corrigés en conséquence, et le processus de propagation et de calcul d'erreur est repris de façon itérative. Le choix d'une fonction d'erreur et la mise en oeuvre de ce type de calcul ne pose pas de difficulté particulière à l'homme du métier.
Un deuxième aspect de la présente invention concerne un procédé d'analyse d'un réseau d'interactions moléculaires dans un système biologique, comportant les étapes suivantes : A') utilisation d'un modèle dynamique du réseau d'interactions moléculaires, ledit modèle étant susceptible d'être obtenu, par un procédé décrit ci-dessus, et construit à partir d'un graphe statique dont les sommets représentent des molécules biologiques du système biologique et les arrêtes représentent des interactions physicochimiques entre ces molécules, et à partir des données expérimentales concernant les taux ou les activités de ces molécules biologiques. C) un état du graphe, mesuré expérimentalement, est choisi comme "état à modifier", et la durée du processus biologique à simuler est définie et découpée en une série de pas de temps, D) plusieurs procédures itératives de simulation sont effectuées, comprenant chacune les étapes suivantes : a) un stimulus est imposé à l'état à modifier, c'est-à-dire que la valeur d'une ou de plusieurs des variables quantitatives associées aux sommets du graphe est modifiée, constituant ainsi un état de départ de la simulation ; b) à partir de l'état de départ de la simulation, un calcul de propagation est effectué au sein du graphe.
Le calcul de propagation au sein du graphe peut être effectué pendant un nombre de pas de temps tel que la durée de la simulation n'excède pas la durée du processus biologique à simuler définie à l'étape C). Toutefois, il est également possible de laisser la simulation se poursuivre au-delà de la durée du processus biologique à simuler définie à l'étape C), par exemple si on cherche à voir si le réseau va à terme trouver un nouvel état stable (état d'équilibre) et si on ne sait pas a priori combien de temps cela va prendre. Il est important de noter que la durée de la simulation définie à l'étape C) peut être plus longue que celle des cinétiques expérimentales utilisées pour le calcul des paramètres (ou plus courte).
Selon une variante du procédé d'analyse d'un réseau d'interactions moléculaires décrit ci-dessus, seules les étapes C), D)a) et D)b) ci-dessus sont effectuées, en utilisant (sans le reconstruire) un modèle dynamique du réseau d'interactions moléculaires choisi, ledit modèle étant susceptible d'être obtenu par un procédé tels que les procédés d'obtention de modèles dynamiques de réseaux d'interactions moléculaires décrits plus hauts.
Un autre aspect particulièrement important de la présente invention est un procédé de sélection de cibles thérapeutiques mettant en œuvre un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, par la mise en œuvre d'un système informatique, comprenant les étapes et caractéristiques suivantes : A') utilisation d'un modèle dynamique du réseau d'interactions moléculaires, ledit modèle étant susceptible d'être obtenu, par un procédé décrit ci-dessus, et construit à partir d'un graphe statique dont les sommets représentent des molécules biologiques du système biologique et les arrêtes représentent des interactions physicochimiques entre ces molécules, et à partir de données expérimentales concernant les taux ou les activités de ces molécules biologiques; C) un état du graphe, mesuré expérimentalement, est choisi comme "état à modifier", et la durée du processus biologique à simuler est définie et découpée en une série de pas de temps, et un état du graphe correspondant à un "état à atteindre" du système biologique est choisi comme "état final du graphe" à atteindre; D) plusieurs procédures itératives de simulation sont effectuées, comprenant chacune les étapes suivantes : a) un stimulus est imposé à l'état à modifier, c'est-à-dire que la valeur d'une ou de plusieurs des variables quantitatives associées aux sommets du graphe est modifiée, constituant ainsi un état de départ de la simulation ; b) à partir de l'état de départ de la simulation, un calcul de propagation est effectué au sein du graphe ; c) un calcul de proximité entre l'"état final du graphe " obtenu à l'issue de l'étape b) et l'état à modifier, ou entre T'état final du graphe " et un état voulu est effectué ; E) à partir de l'ensemble des proximités statistiques calculées à l'étape D), les sommets, et les stimuli imposés sur ces sommets, sont hiérarchisés, les sommets hiérarchisés correspondant à des cibles thérapeutiques classées.
Bien entendu, et comme pour le procédé d'analyse d'un réseau d'interactions, le procédé de sélection de cibles thérapeutiques selon l'invention peut être mis en œuvre en effectuant uniquement les étapes C) à E) ci-dessus en utilisant, sans le reconstruire, un modèle dynamique susceptible d'être obtenu par des méthodes d'obtention de tels modèles, décrites plus haut. De même, l'étape D) b) peut être poursuivie au-delà de la durée spécifiée à l'étape C)
L'étape A') des procédés ci-dessus peut être réalisée de la même façon que les étapes A) et B) des procédés d'obtention de modèles dynamiques de réseaux d'interactions décrits plus haut.
Dans ces procédés, l'étape C) peut être réalisée en tenant compte des éléments suivants, lorsque le cas de figure s'y prête : Pas de temps : La durée du processus biologique à simuler est découpée en une série de pas de temps, régulièrement espacés ou non ; les pas de temps sont définis de façon à être préférablement plus petits que les durées expérimentales réelles séparant les séries de données expérimentales quantitatives utilisées pour le calcul des paramètres des relations. La définition de ces pas de temps est rendue nécessaire par le fait que tout processus informatique de simulation dynamique consiste à calculer des états à des temps discrets, rendant la discrétisation du temps nécessaire. On obtient donc une série de temps consécutifs, sur lesquels va être effectuée la simulation. Le premier temps de la série est appelé le temps initial. Ce temps initial correspond à l'état de départ du graphe, défini plus bas.
Etats du graphe pour les simulations : Un état du graphe, mesuré expérimentalement, et correspondant à un état que l'on veut modifier du système biologique, est défini (par exemple un état pathologique). Cet état est appelé "l'état à modifier". Dans certains cas, l'homme du métier peut savoir que les différences entre l'état à modifier et l'état stable de référence concernent essentiellement un sous-ensemble des molécules du réseau, et décider de ne mesurer expérimentalement que les valeurs des variables correspondantes, les autres Xj étant alors, par défaut, fixés aux valeurs de l'état stable de référence. Un état du graphe (mesuré expérimentalement ou défini arbitrairement), correspondant à un état que l'on veut atteindre du système biologique, est éventuellement défini (par exemple un état sain). Cet état est appelé "l'état à atteindre".
Identification de molécules-cibles thérapeutiques pour une pathologie donnée :
Dans une mise en œuvre préférentielle de l'invention, l'état à modifier et l'état à atteindre sont définis comme suit : On pratique les simulations d'actions sur les sommets du graphe à partir d'un état à modifier du graphe identique ou similaire à son état tel qu'observé expérimentalement dans la condition pathologique (par exemple par criblage d'expression des ARN messagers sur puces à ADN à partir de tissus pathologiques).
On définit l'état à atteindre comme étant un état proche d'un état non pathologique de référence (tel que mesuré lui aussi par l'observation expérimentale de la condition non pathologique, par exemple par criblage d'expression des ARN messagers sur puces à ADN à partir de tissus sains).
Le processus de simulation consiste alors à identifier les sommets, et les stimuli sur ces sommets, qui, en partant de l'état à modifier (l'état pathologique), permettent le mieux de faire évoluer le graphe (en partie ou entièrement) vers un état proche de l'état à atteindre (état non pathologique).
Identification des molécules-cibles thérapeutiques de traitements existants ou en cours de développement, et pour lesquels aucune ou une partie seulement des cibles sont connues (ce qui est le cas de nombreux médicaments actuels).
Dans ce cas, l'état à modifier est défini comme ci-dessus, et l'état à atteindre est défini comme l'état, ou un état proche, de celui obtenu expérimentalement lors de l'administration de ce traitement (tel que mesuré par exemple par criblage d'expression des ARN messagers sur puces à ADN à partir de tissus pathologiques qui ont été soumis au traitement concerné).
Le processus de simulation consiste alors à identifier les sommets, et les stimuli sur ces sommets, qui, en partant de l'état à modifier (l'état pathologique), permettent le mieux de faire évoluer le graphe (en partie ou entièrement) vers un état proche de l'état à atteindre (état pathologique sous traitement). Cette mise en œuvre particulière peut aussi être réalisée en définissant l'état à modifier comme tout état G possible du système biologique étudié (par exemple l'état sain), et l'état à atteindre comme l'état obtenu après l'administration du traitement concerné au système biologique à l'état G.
Dans les procédés d'analyse de réseaux d'interactions et de sélection de cibles selon l'invention, l'étape D) est réalisée en considérant les éléments suivants :
Stimulus : Un stimulus est imposé à l'état à modifier. Ce stimulus est exercé sous la forme de la variation de la valeur d'une ou de plusieurs variable(s) du graphe (correspondant à un ou plusieurs sommet(s)), c'est-à- dire d'une augmentation ou d'une diminution de cette ou ces valeur(s), selon la simulation souhaitée. Les valeurs de toutes les autres variables restent inchangées. On obtient donc un nouvel état du graphe, qui est
"l'état de départ' de la simulation. L'état de départ et l'état à modifier ne diffèrent donc que par la valeur de la ou des variable(s) modifiée(s), toutes les valeurs de toutes les autres variables étant identiques. Cet état est défini comme correspondant au premier temps de la simulation. Dans une mise en œuvre particulière du procédé, les stimuli portent à chaque fois sur un seul sommet.
Propagation : A partir de l'état de départ de la simulation, un calcul de propagation est effectué au sein du réseau. Cette propagation consiste à calculer les nouvelles valeurs de toutes les variables au pas de temps suivant, aboutissant à un nouvel état du graphe, et à recommencer le calcul à partir de ce nouvel état pour le pas de temps suivant, et ainsi de suite. Cette propagation se prolonge pendant le nombre de pas de temps (donc la durée biologique) définie par l'expérimentateur en fonction de la question biologique posée. Elle peut éventuellement être prolongée jusqu'à l'apparition d'un nouvel état stable du graphe (un nouvel état d'équilibre), ou être arrêtée avant. Au terme de cette simulation, un nouvel état ("état final') du graphe est obtenu.
Itération : Le processus précédent est répété avec un nouveau stimulus, portant sur un ou plusieurs autre(s) sommet(s) du graphe, ou portant éventuellement sur le(s) même(s) sommet(s) du graphe avec l'imposition d'une nouvelle valeur à la ou aux variable(s).
Ce processus peut être répété de façon itérative sur tous les sommets individuellement, éventuellement en imposant plusieurs valeurs (en nombre fini) par variable de manière à tester des gammes d'activation ou d'inhibition sur tous les nœuds. Dans ce cas, le résultat de l'étape E) est une hiérarchisation des sommets, et des stimuli imposés sur ces sommets. Ce classement correspond donc au classement des sommets, de celui sur lequel un stimulus est le plus susceptible d'aboutir à l'état voulu à partir de l'état à modifier, jusqu'à celui sur lequel un stimulus est le moins susceptible d'avoir cet effet. A chaque proximité correspond en effet un et un seul sommet et une et une seule valeur de stimulation sur ce sommet.Si l'effet recherché est l'amélioration d'un état pathologique, ce classement est celui des cibles thérapeutiques potentielles, de la plus probable à la moins probable.
Bien que présenté ici de manière séquentielle, l'ensemble des propagations effectuées peut être calculé de manière parallèle.
A l'étape D)c), la proximité de chaque état final obtenu à l'étape D)b) peut être calculée soit par rapport à l'état à modifier choisi à l'étape C), soit par rapport à un autre état, mesuré expérimentalement ou déterminé arbitrairement, et défini comme P"état à atteindre", qui peut être, par exemple, un état sain. Il peu s'agir de l'état de référence défini plus haut.
Une fois les calculs de proximité de graphes effectués pour toutes les simulations, l'étape E) consiste à classer l'ensemble des multiplets (sommet(s) du graphe - stimulus) en ordre hiérarchique (croissant ou décroissant) correspondant directement à l'ordre hiérarchique (croissant ou décroissant, respectivement) des proximités qui leur sont associées. Aux sommets du graphe correspondent directement les molécules du réseau biologique, qui sont donc hiérarchisées de fait.
Cette hiérarchisation ne pose aucun problème technique à l'homme de l'art, les proximités étant des valeurs numériques positives pouvant être directement hiérarchisées de la plus grande à la plus petite, ou inversement. Le résultat de ce classement peut être avantageusement produit sous forme de liste ou de tableau, ou sous tout autre type de format, et / ou stocké dans une base de données en vue d'une utilisation ultérieure.
Quels que soient les niveaux de proximité des graphes, une hiérarchisation des multiplets (sommet(s) du graphe - stimulus) selon cette méthode sera toujours obtenue. L'invention permet donc toujours d'obtenir un résultat, en fonction des connaissances biologiques et des techniques de mesure expérimentales utilisées. Elle ne requiert pas de connaissance préalable étendue des processus physiopathologiques moléculaires en œuvre dans le processus pathologique analysé. Toutes les molécules du réseau d'interactions moléculaires sont considérées a priori (avant mise en œuvre de l'invention) comme des molécules cibles thérapeutiques potentielles sans en exclure aucune, les molécules cibles thérapeutiques étant sélectionnées a posteriori (après mise en œuvre de l'invention) sur des critères statistiques objectifs (calculs de proximités). Cette méthode est utilisable de façon systématique et automatisée quelle que soit la pathologie étudiée, dès lors qu'il est possible de définir un état à modifier. Ceci la rend notamment particulièrement adaptée à une utilisation dans le cadre de processus industriels de sélection systématique à grande échelle de cibles thérapeutiques, en utilisant les données expérimentales fournies par les technologies de criblages moléculaires à grande échelle. Dans le cas de l'identification de cibles thérapeutiques, le classement hiérarchique des molécules du réseau biologique correspond directement au classement hiérarchique de ces molécules considérées comme cibles thérapeutiques. L'invention permet donc d'obtenir un classement des cibles thérapeutiques potentielles hiérarchisées selon des critères statistiques objectifs, en fonction des données expérimentales (mesures des Xj) et des connaissances fonctionnelles du réseau (existence d'interactions moléculaires). Dans les cas où il est possible de définir à la fois un état à modifier et un état à atteindre, les meilleures cibles thérapeutiques potentielles sont considérées comme étant celles correspondant aux proximités les meilleures avec l'état à atteindre.
Dans les cas où la définition d'un état à atteindre n'est pas possible (ce qui devrait être exceptionnel, l'état sain pouvant à priori toujours être utilisé par défaut comme état à atteindre pour les processus pathologiques), il est possible de hiérarchiser les multiplets (sommet(s) du graphe - stimulus) par leur proximité avec l'état à modifier, et de classer les molécules du réseau biologique considérées comme cibles thérapeutiques potentielles en suivant une hiérarchie directement inverse de celle des proximités : les meilleures cibles thérapeutiques potentielles sont considérées comme étant celles correspondant aux proximités les plus mauvaises par rapport à l'état à modifier.
Un point important est que cette invention permet non seulement d'identifier des molécules-cibles thérapeutiques, mais aussi de prédire le sens de l'action thérapeutique qu'il sera nécessaire d'appliquer sur ces molécules (activation ou inhibition).
Les cibles thérapeutiques sont donc sélectionnées à partir des données concernant l'ensemble des molécules étudiées, et non seulement celles concernant spécifiquement les molécules-cibles, puisque le critère utilisé pour la hiérarchisation dépend de l'évolution du graphe dans son ensemble, donc de l'ensemble des mesures expérimentales d'expression et/ou d'activation de toutes les molécules représentées dans le graphe, et non la simple évolution des mesures expérimentales d'expression et/ou d'activation des seules molécules cibles. Il s'agit donc bien d'une méthode integrative répondant aux besoins actuels tels que définis plus haut, notamment en ce qui concerne des maladies à déterminisme multi-factoriel, apportant clairement un progrès par rapport aux méthodes de sélection de cibles thérapeutiques existantes.
La méthode d'identification des cibles décrite ci-dessus comporte les caractéristiques avantageuses suivantes : - Les calculs sont fondés sur méthode non probabiliste, ce qui élimine toute limitation en termes de temps de calcul, au contraire des méthodes des équations stochastiques et des réseaux bayésiens. - L'invention intègre les données expérimentales quantitatives, ce qui la différencie des méthodes qualitatives (réseaux booléens, formalismes logiques généralisés, formalismes fondés sur des règles), permet d'éviter des contraintes et hypothèses sur le fonctionnement du réseau, et permet d'augmenter la fiabilité des simulations. - Le fait de définir les variables comme similaires aux données expérimentales effectivement mesurables permet de calculer les paramètres des relations de façon optimale (sans avoir à extrapoler les valeurs des variables). - Le fait d'établir, pour tout sommet du graphe, une relation directe entre la variable qui lui est associée et les variables associées aux sommets du graphe agissant sur ce sommet permet la mise en œuvre directe de méthodes de calcul des paramètres dérivées des méthodes d'apprentissage de réseaux de neurones par calcul de l'erreur minimale qui sont compatibles avec des réseaux de grande taille en termes de temps de calcul. - Une fois les paramètres calculés, les simulations sont très peu coûteuses en temps de calcul, le réseau étant déterministe. Ceci est aussi compatible avec l'application de l'invention à des réseaux de grande taille. - Les limitations de divergence introduites dans les relations ou fonctions permettent de pratiquer des simulations sur des cinétiques longues et des réseaux de grande taille avec une fiabilité satisfaisante. - Les connaissances de l'existence d'interactions entre les molécules du réseau, et de l'orientation d'une partie de ces interactions, sont suffisantes pour la mise en œuvre de l'invention. La connaissance du type d'interaction (activation ou inhibition) peut être avantageusement utilisée lorsqu'elle est disponible, mais elle n'est pas indispensable. Aucune autre connaissance qualitative supplémentaire concernant le réseau n'est requise. Pour les grands réseaux d'interactions moléculaires (plus d'une centaine de molécules) ces connaissances sont en général les seules disponibles aujourd'hui. La qualification de cette méthode, suivant les critères considérés dans les tableaux 1 et 2 ci-dessus, est donc la suivante :
Figure imgf000048_0001
Tableau 3 II est important de noter que le formalisme permettant de prendre en compte l'inertie/ tendance au retour à l'état initial est spécifique à l'invention. En effet, dans la méthode de la présente invention, les conséquences des interactions sont représentées comme résultant d'une résistance au changement des taux de molécules suite à une modification quantitative de l'activité biologique d'au moins une molécule interagissant sur elles et une tendance à revenir à l'état initial ; cette représentation permet d'éviter de faire des hypothèses sur le fonctionnement du système (effets de seuil, types de réactions chimiques, etc.) et de tenir compte des données ou variables éventuellement non connues ou non mesurées, l'inertie et la tendance au retour à l'état initial représentant la résultante des multiples phénomènes biologiques impliqués dans une interaction donnée (temps de synthèse de la molécule, existence d'un rétro-contrôle négatif concomitant, temps de transport des molécules jusqu'au compartiment cellulaire où elles sont actives, etc.) ; le formalisme de l'invention est donc fondamentalement différent de celui des autres méthodes existantes (comparer avec le tableau 1).
Figure imgf000049_0001
Tableau 4
Selon une variante des procédés de sélection de cibles thérapeutiques décrits plus haut, un premier classement hiérarchique des sommets, considérés individuellement, est obtenu en effectuant les étapes A) à E) en imposant, pour chacune des simulations de l'étape D), des stimuli qui concernent un sommet unique ; une étape D2) est ensuite effectuée, correspondant à l'étape D) dans laquelle les stimuli imposés à chaque simulation sont exercés sur deux sommets, soit en testant toutes les combinaisons de deux sommets possibles, soit en limitant ces calculs aux combinaisons de deux sommets parmi un certain nombre des sommets les mieux classés à l'étape E). Enfin, une étape E2) de classement hiérarchique des associations de deux sommets sur lesquels des stimuli sont le plus susceptibles d'avoir l'effet voulu est effectuée à partir de l'ensemble des proximités statistiques calculées à l'étape D2).
A partir de la variante ci-dessus, les étapes D) et E) peuvent être répétées de façon itérative, en augmentant à chaque fois le nombre de sommets sur lesquels sont exercés les stimuli. Ainsi, le procédé peut comporter une étape D3) suivant l'étape E2) et correspondant à l'étape D) dans laquelle les stimuli imposés à chaque simulation sont exercés sur trois sommets, soit en testant toutes les combinaisons de trois sommets possibles, soit en limitant ces calculs aux combinaisons de trois sommets choisis parmi un certain nombre des sommets les mieux classés à l'étape E) et des combinaisons de deux sommets les mieux classées à l'étape E2), ladite étape D3) étant suivie d'une étape E3) de classement hiérarchique des associations de trois sommets sur lesquels des stimuli sont le plus susceptibles d'avoir l'effet voulu. Des étapes D4) et E4), avec des stimuli sur 4 sommets peuvent ensuite être rajoutées, et ainsi de suite. Ces procédés de sélection de cibles thérapeutiques comportent de préférence une étape finale de classement statistique des proximités de graphes de toutes les simulations effectuées, intégrant l'ensemble des classements précédemment obtenus.
Dans les procédés de l'invention, lorsqu'une simulation implique des stimuli sur une combinaison de sommets, les stimuli exercés sur ces différents sommets peuvent être appliqués simultanément ou non, c'est-à-dire que la simulation peut tenir compte d'un décalage temporel entre les stimuli exercés sur différents sommets. Dans une mise en œuvre de l'invention, la relation entre Xj et les Xj est établie, pour au moins une partie des interactions physico-chimiques entre les molécules du réseau, par une relation inertielle découlant de celle de l'oscillateur harmonique en physique, associée à un facteur d'amortissement suffisamment important pour limiter l'oscillation à un seul cycle.
Plus précisément, une relation de ce type entre Xj et chaque Xj, deux à deux est de la forme : Wjj .Xj = mi .(d2Xi / dt2) + 2 .λ .(dXi / dt) + ωy 2.Xj, dans laquelle mi .(d2Xi / dt2) + Guy2.Xj correspond au terme inertiel (i), 2 .λ .(dXj / dt) correspond au terme de retour à l'état initial (ii), Xi est une variable associée à la molécule i dXi / dt est la dérivée de Xj en fonction du temps d2X, / dt2 est la dérivée seconde de Xj en fonction du temps Xj est une variable associée à la molécule j, ιτii représente l'inertie de i, λy régit le retour à l'état d'équilibre de Xj, la pulsation ωy correspond au temps de réponse de Xj à la variation de Xj, et Wjj est un facteur de couplage représentant la force de l'interaction entre les molécules i et j, correspondant à une pondération de l'effet de chaque molécule j sur la molécule i vis-à-vis de la résultante de l'ensemble des effets combinés de toutes les molécules j exerçant un effet sur i.
Selon une autre mise en œuvre des procédés de l'invention, pour au moins une partie des interactions physico-chimiques entre les molécules du réseau, la relation entre les variables Xj et Xj, deux à deux est établie par une relation sigmoïde comportant un facteur de retardement associée à une fonction de décroissance linéaire. Un autre type de relation entre les variables Xj et Xj, décrit plus en détail ci- après, utilisable dans les procédés de l'invention pour modéliser au moins une partie des interactions physico-chimiques entre les molécules du système biologique, est de la forme : (dXj/dt) = Kn . [ 1 / (1 + e ∑ wij Xj -bi) ] - K2i . Xj ; où : le terme sigmoïde K-π . [1 / (1 + e"∑ lj Xj " bi)] correspond au terme inertiel (i), et le terme K2j . Xj correspond au terme de retour à l'état initial (ii), avec : Xj = variable associée au sommet i, Xj = variable associée au sommet j, wij = facteur de couplage représentant la force de l'interaction entre les molécules i et j, correspondant à une pondération de l'effet de chaque molécule j sur la molécule i vis-à-vis de la résultante de l'ensemble des effets combinés de toutes les molécules j exerçant un effet sur i., bi = facteur de retardement, Kn = facteur de limite maximale de variation de Xi, et K2j = facteur de retour à l'équilibre.
Dans les procédés ci-dessus, la relation entre les variables Xi et Xj, peut également être, pour au moins une partie des interactions considérées, une fonction polynôme de type wy Xj = Σ bmj.Xim = typ-Di .Xip"1 + ... + b3, .Xi3 + b2i .Xj2 + bii .Xi + boi ,
d'ordre strictement inférieur au nombre p de couples (Xit, Xjt) de valeurs expérimentales du niveau de taux ou d'activité Xj ou Xj des molécules i et j, respectivement, à différents instants t, les paramètres bmj étant calculés à partir des p couples expérimentaux (Xit, Xjt) disponibles, et w étant un facteur de couplage représentant la force de l'interaction entre les molécules i et j, correspondant à une pondération de l'effet de chaque molécule j sur la molécule i vis-à-vis de la résultante de l'ensemble des effets combinés de toutes les molécules j exerçant un effet sur i.
Des fonctions de type dérivée :
Figure imgf000053_0001
[m-0-->p'-1] p' étant un entier tel que 1 < p' > p - 1 , et p étant défini tel que ci-dessus, peuvent également être utilisées dans les procédés de l'invention pour modéliser au moins une partie des interactions physico-chimiques entre les molécules du système biologique.
Ceci peut notamment être mis en œuvre avec p'=3.
La résultante globale de n interactions exercées par des molécules 1 à n sur une molécule i peut être, dans les procédés de l'invention, et pour au moins une partie des molécules du réseau, une somme pondérée des actions des molécules 1 à π sur la molécule i, de la forme FG (∑ j -» i ) = ∑ a . i, où [j.1- n] 0:1-->n] fji est la fonction associée à l'arc (i, j) pour chaque couple (i, j) et ay = (dXj/dt) / Σ (dXj/dt). B-1-»n]
Une telle somme pondérée peut également être faite avec ay = (d2Xj/dt2) (d2Xj/dt2).
Figure imgf000053_0002
La présente invention porte également sur un procédé de détermination du mode d'action d'un xénobiotique, consistant à mettre en œuvre un procédé d'analyse d'un réseau d'interactions moléculaires dans un système biologique, tels que ceux décrits plus haut, dans les conditions suivantes : (i) le système biologique dans lequel un réseau d'interactions moléculaires est étudié est concerné par l'action du xénobiotique ; (ii) H'état à modifier" choisi à l'étape C), correspond à un état observé expérimentalement avant l'administration dudit xénobiotique ; (iii) on identifie les modifications à apporter au cours de l'étape D)a) pour que le calcul effectué à l'étape D)b) montre une évolution du système vers un état proche de l'état observé après administration du xénobiotique.
Un autre aspect de l'invention est une méthode de prédiction d'effets indésirables de traitements appliquant un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, par la mise en œuvre d'un système informatique.
Dans cet aspect de l'invention, les étapes et caractéristiques de la méthode sont les mêmes que précédemment, la seule modification consistant dans l'adaptation suivante :
Une fois identifiées les molécules-cibles d'un traitement, on analyse sur des parties du graphe représentatives de fonctions physiologiques connues, par des simulations mettant en jeu la même méthode que dans les aspects précédents de l'invention (étapes A à E, éventuellement A à Ek lorsque les étapes D et E sont répétées de façon itérative en appliquant les stimuli sur des combinaisons de sommets allant jusqu'à k sommets), les conséquences de l'application du traitement sur ces molécules cibles. Cette analyse consiste à identifier les éventuelles évolutions de ces sous-parties de graphes vers de nouveaux états proches d'autres états pathologiques de référence (tels que définis par l'observation expérimentale de ces conditions pathologiques, selon des méthodes similaires à ce qui est décrit plus haut).
A titre d'exemple, l'observation lors des simulations de l'évolution du sous graphe de l'apoptose vers un état final ayant une grande proximité avec un état de référence de ce graphe correspondant à une activation de cette voie physiologique (telle que définie à partir de données concernant un ou des tissus affectés par des processus de dégénérescence cellulaire) permet de prédire un effet de toxicité cellulaire du traitement dans le ou les tissu concernés.
Cet aspect de l'invention consiste donc à mettre en œuvre un procédé d'analyse tel que décrit plus haut, dans les conditions suivantes : (i) le système biologique dans lequel un réseau d'interactions moléculaires est étudié est concerné par le traitement ; (ii) les modifications de l'étape D)a) correspondent aux modifications des niveaux de taux ou d'activité des molécules cibles observées ou souhaitées lors de l'application du traitement ; (iii) l'étape D)b) de calcul de l'évolution du système biologique est suivie d'une analyse de sous-parties du système correspondant à des fonctions physiologiques connues, afin d'identifier les éventuelles évolutions de ces sous-parties vers des états proches d'états pathologiques de référence.
La présente invention porte également sur un procédé pour hiérarchiser des cibles thérapeutiques potentielles pour une pathologie, consistant à identifier des cibles thérapeutiques par un procédé selon l'invention, puis à prédire les éventuels effets indésirables d'un traitement visant ces cibles, et enfin à déterminer le rapport "bénéfice thérapeutique / effets indésirables" d'une action sur chacune des cibles thérapeutiques potentielles.
Comme exposé ci-dessus, un des avantages principaux de la présente invention, dans ses différents aspects, est de permettre de travailler sur des graphes ou réseaux de molécules en interaction comportant un grand nombre de molécules. Dans l'ensemble des procédés de l'invention, décrits plus haut, le nombre de variables Xi du réseau d'interactions moléculaires considéré est donc de préférence supérieur à 100, supérieur à 200, voire supérieur à 300. L'invention concerne aussi un procédé d'analyse tel que décrit plus haut faisant appel à l'utilisation des réseaux d'interaction moléculaire de l'invention, lesdits réseaux étant associés pour former un hypergraphe de réseaux.
Selon cette variante de réalisation de l'invention, le nombre de variables Xi de chaque réseau d'interactions moléculaires est inférieur à environ 100 et le nombre de réseaux associés pour former l'hypergraphe est compris entre 2 et environ 100.
Un autre aspect de l'invention est une méthode d'extension des graphes à partir de résultats de criblages expérimentaux des variations de taux d'expression ou d'activité de molécules, appliquant un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, par la mise en œuvre d'un système informatique.
Dans cet aspect de l'invention, les étapes et caractéristiques de la méthode sont les mêmes que précédemment, la seule modification consistant dans l'adaptation suivante : Dans cette application, la méthode est mise en oeuvre pour identifier de nouvelles interactions moléculaires. Ceci peut être réalisé par le couplage de la méthode de l'invention décrite plus haut, avec des méthodes statistiques de recherche de corrélation entre des points dans un espace à n dimensions (par exemple analyse factorielle, classifications hiérarchiques, etc.) telles que (mais de façon non exclusive) celles utilisées à ce jour pour rechercher des corrélations de l'expression de gènes à partir des résultats de criblage d'ARN messagers sur puces à ADN ("clustering" de gènes). A titre d'exemple de méthode de "clustering", on peut citer Eisen MB, Spellman PT, Brown PO and Botstein D (1998), Cluster Analysis and Display of Genome-Wide Expression Pattems, Proc Natl Acad Sci U S A 95, 14863-8. Un exemple de système logiciel d'accès libre permettant de réaliser des analyses de clustering disponible sur internet est le logiciel "cluster 3.0", développé par le Laboratory of DNA Information Analysis of Human Génome Center, http://www.ims.u-tokyo.ac.ip/imswww/index- e.htmllnstitute of Médical Science, Universitv of Tokyo, au Japon (4-6-1 Shirokanedai, Minato-ku, Tokyo 108-8639 JAPAN). Le logiciel "cluster 3.0" est disponible sur le site internet http://bonsai.ims.u- tokyo.ac.jp/~mdehoon/software/cluster/. Les données expérimentales utilisées peuvent par exemple être celles produites par les criblages d'expression d'ARN messagers sur puces à ADN. Ce couplage consiste à utiliser les paramétrages calculés par la mise en œuvre de l'invention pour re-calculer une nouvelle matrice de données expérimentales de mesure d'expression de taux ou d'activité des molécules, en éliminant des matrices de résultats expérimentaux d'origine les facteurs d'interactions moléculaires inclus dans le modèle dynamique paramétré (tels que la composante de résistance dynamique ou inertielle), puis à effectuer les recherches de corrélation. Ce "nettoyage" des matrices de résultats d'origine consiste en d'autres termes à en éliminer le "bruit statistique" lié à ces facteurs, ces facteurs étant alors considérés comme introduisant des distorsions, dans les mesures réellement observées des taux d'expression ou d'activité des molécules, par rapport à ce qu'auraient été ces mesures, d'un point de vue théorique, en l'absence de ces facteurs.
A titre d'exemple, la résistance dynamique de l'expression d'un gène A donné (donc l'inertie de la modification du taux d'ARN messager correspondant) à deux stimulations distinctes exercées par les molécules B et C (elles mêmes distinctes) peut varier, empêchant avant tout "nettoyage" de ce type de mettre en évidence à la fois une corrélation entre l'expression de A et l'activité de la molécule B, et une corrélation entre l'expression de A et l'activité de la molécule C.
L'invention porte donc sur l'utilisation d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique susceptible d'être obtenu par un procédé tel que décrit plus haut, pour étendre un graphe statique dont les sommets représentent des molécules biologiques et les arcs représentent des interactions physico-chimiques entre ces molécules.
D'autres avantages et caractéristiques de la présente invention apparaissent dans les exemples ci-après de mise en œuvre pratique des procédés de l'invention, qui illustrent de façon non limitative les méthodes décrites ci-dessus.
Les schémas et figures ci-après illustrent également certains aspects de l'invention :
La figure 1 représente un schéma synthétique des diverses étapes d'une méthode d'identification des cibles selon la présente invention. La figure 2 représente un schéma du graphe construit dans l'exemple 4. Ce graphe comporte 116 molécules dans le réseau d'interactions moléculaires (116 sommets), et 329 interactions moléculaires entre ces molécules. Chaque rectangle représente une molécule du réseau (= un sommet du graphe). Chaque flèche représente une interaction entre deux molécules (= une arrête du graphe) ; le sens des flèches représente le sens des interactions : si la flèche va de la molécule A vers la molécule B, cela signifie que la molécule A a une action d'activation ou d'inhibition potentielles de la molécule B ; certaines flèches sont à double sens : interaction bilatérale. Le texte au sein de chaque rectangle correspond à l'abréviation du nom de la protéine telle que décrite dans le texte. Le calcul des paramètres des relations entre les sommets du graphe et les simulations ont été réalisées sur l'ensemble de ce graphe (exemple 4).
La figure 3 montre des exemples graphiques de cinétiques calculées (triangles) et observées (carrés), pour quelques gènes (exemple 4). Figure 3A : ORF YBL015W (ACM). Figure 3B : ORF YMR169C (ALD3). Figure 3C : ORF YIL125W (KGD1). Figure 3D : ORF YNL071W (PDA2). Figure 3E : ORF YAL054C (ACS1). Figure 3F : ORF YFL018C (LPD1).
La figure 4 représente un schéma du graphe construit dans l'exemple 5. Ce graphe comporte 133 molécules dans le réseau d'interactions moléculaires (133 sommets), et 407 interactions moléculaires entre ces molécules.
La signalétique des rectangles, flèches et des textes au sein de chaque rectangle est la même que celle précédemment décrite, en référence à la figure 2.
La figure 5 montre des exemples de courbes de paramétrage, dans lesquels les cinétiques mesurées expérimentalement sont représentées en blanc et les cinétiques calculées par simulation sont représentées en noir, pour quelques molécules (exemple 5). Figure 5A : ICL 1 (YER065C).
Figure 5B : IDH1 (YNL037C). Figure 5C : ACM (YBL015W). Figure 5D : PCK1 (YKR097W).
La figure 6 représente un schéma de classification des molécules du réseau par classification hiérarchique des distances calculées entre d'une part l'état à atteindre et d'autre part les états obtenus par simulation
(exemple 5).
Les ordonnées correspondent aux valeurs de distance calculées. En abscisse les 133 molécules du réseau sont classées de gauche à droite de celle associée à la distance la plus faible à celle associée à la distance la plus élevée, chaque point correspondant à une molécule du réseau.
Exemples
Exemple 1 : Mise en œuyre pratique de l'étape A)
- (1) Les variables associées aux sommets du graphe
Soit i une molécule donnée du réseau, et Xj sa quantité (ou sa concentration) au sein du système biologique étudié. Soit xi0 la mesure expérimentale effectivement réalisée de i à un « état étalon » du système biologique, utilisé lors des mesures. Soit XH la mesure expérimentale effectivement réalisée de i à un instant t. La variable utilisée est :
Figure imgf000060_0001
L'éfaf étalon est un état mesurable utilisé pour pratiquer les mesures biologiques, contre lequel toutes les autres mesures sont quantifiées. Il peut correspondre à un état artificiel du système, par exemple à un pool de plusieurs échantillons biologiques prélevés dans des conditions expérimentales différentes (état artificiel), ou à un état réellement observable (non artificiel) du système. Cette variable correspond bien au type de mesures biologiques effectivement réalisables. A. titre d'exemple, lors des mesures de taux d'ARN messagers sur puces à ADN, la mesure effectivement réalisée pour chaque ARN à un temps expérimental t donné est le rapport du signal émis par l'hybridation dès ARN présents dans l'échantillon biologique au temps t sur le signal émis par les ARN de même type présents dans l'échantillon à un état étalon du système biologique étudié (par exemple le temps initial de l'expérience biologique). Seule cette mesure peut être considérée comme fiable, la quantité réelle de molécules d'ARN n'étant pas directement mesurable car elle dépend de paramètres expérimentaux non directement contrôlés (rendement des réactions de marquage des sondes, rendement des hybridations sur la puces, etc., ces paramètres différant de façon non prédictible entre deux ARN de type différent donnés). La quantité de signal mesuré à l'état étalon sert donc d'étalon de mesure pour celle aux autres temps, en se fondant sur l'hypothèse que pour un type d'ARN donné, les paramètres expérimentaux influant sur le signal finalement émis sont les mêmes.
Xi correspond donc directement aux mesures quantitatives biologiques réellement productibles dans l'état actuel des techniques de biologie moléculaire. Les variables Xj, Xj etc. sont donc égales à (xit/xio), ( jt/ jo) etc. - (2) Les relations associées aux arrêtes du graphe et reliant les variables :
Soient n sommets ji, j2,..., jn du graphe (correspondant à n molécules du réseau) qui agissent sur un sommet i (orientation du graphe des j vers i).
Ces relations définissent une relation directe entre Xj et les Xj (Xj-i,
Figure imgf000061_0001
Terme inertiel de ces relations : Ce terme correspond à une fonction continue des Xj. Ce terme comporte une composante inertielle. Par inertie, on entend le fait que Xj présente une résistance au changement suite à une variation des Xj : plus précisément, ce terme de la relation doit permettre de rendre compte du comportement suivant des variables : suite à une variation donnée d'un ou plusieurs des Xj, la vitesse de variation de Xi va être initialement faible, puis s'accélérer progressivement.
Ce terme doit aussi permettre de rendre compte du comportement suivant des variables : suite à la variation d'un ou plusieurs des Xj, Xj va progressivement atteindre une nouvelle valeur finie correspondant à la variation maximale de Xj (pic de variation) ; ceci revient à dire que la vitesse de variation de Xj, après avoir augmenté, va diminuer et progressivement tendre vers 0. Il y a donc une inflexion de la courbe de Xj en fonction du temps. Commentaires :
Le fait de comporter une composante inertielle introduit de fait l'expression d'un retard temporel de la variation de Xi suite à la variation de Xj : en l'absence d'autres interactions s'exerçant sur i, le pic de variation de Xj tend à survenir après le pic de variation de Xj. Le fait de comporter une composante inertielle permet donc de rendre compte du décalage temporel dans les variations des Xj lors de la propagation des activation / inhibitions dans le réseau. A l'inverse, le fait d'introduire un simple décalage temporel par d'autres méthodes mathématiques n'introduira pas systématiquement un terme inertiel.
Terme de retour à l'état initial de ces relations : Ce terme tend à ramener Xj à son niveau initial.
Il correspond à une fonction continue de Xj toujours croissante en valeur absolue avec Xj : plus la variation de Xj est forte, plus ce terme augmente en valeur absolue et tend à ramener Xj à son niveau de départ.
Le couplage de ces deux termes permet de définir une relation générale pouvant rendre compte de cinétiques variables et non linéaires.
Plusieurs relations mathématiques peuvent comporter ces caractéristiques et être utilisées pour établir une relation entre Xj et les Xj. Le fait de présenter une inflexion, d'être contraint par une limite maximale finie, et de tendre au retour à l'état initial peut être obtenu notamment par :
Dans un aspect de l'invention, la relation entre Xj et les Xj est établie par une relation inertielle découlant de celle de l'oscillateur harmonique en physique associée à un facteur d'amortissement suffisamment important pour limiter l'oscillation à un seul cycle.
Pour plus de clarté dans la description, une relation de ce type entre Xj et chaque Xj, deux à deux est de la forme : (1 ) wy .Xj = mi .(d2Xi / dt2) + 2 .λ0 .(dXj / dt) + ωy2 .Xi Le terme : rrii .(d2Xi / dt2) + ωy2 .Xj correspond au terme inertiel Le terme : 2 .λy .(dXj / dt) correspond au terme de retour à l'état initial
Xi = variable associée à la molécule i dX| / dt = dérivée de Xj en fonction du temps d2Xj / dt2 = dérivée seconde de Xj en fonction du temps
Xj = variable associée à la molécule j mj = inertie de i (résistance au changement) λy = amortissement (régit le retour à l'état d'équilibre de Xj ). ω = pulsation (temps de réponse de Xj au stimulus représenté par Xj) wy = facteur de couplage représentant la force de l'interaction entre les molécules i et j, correspondant à une pondération de l'effet de chaque molécule j sur la molécule i vis à vis de la résultante de l'ensemble des effets combinés de toutes les molécules j exerçant un effet sur i.
Une variante formellement équivalente également utilisable de cette relation entre Xj et Xj consiste à se ramener à un cas. particulier où mi est considéré comme une constante ; dans ce cas on peut éliminer le paramètre mj de l'équation. La formulation de l'équation reste globalement la même : wy .Xj = (d2Xj / dt2) + 2 .λjj (dXj / dt) + ωy2.Xi
Dans les deux cas, si : | λy | > 1 , l'amortissement est tel qu'il n'y a plus qu'une seule « oscillation » de Xj. En d'autres termes, cette relation fait varier Xi à partir de son niveau de départ, jusqu'à une variation maximale, puis tend à le faire revenir à son état initial.
Pour mettre en œuvre l'invention, on définit donc une relation entre Xj et l'ensemble des Xj :
Dans un aspect de l'invention, celle-ci est définie par une sommation pondérée des effets des Xj sur la résultante sur Xj : (2) ∑ (wy . Xj) + q = mi .(d2Xi / dt2) + 2 .λy .(dXi / dt) + ωB 2 .Xi D-1 n]
Ou :
(3) ∑ (wy . Xj) + Q =(d2Xi / dt2) + 2 .Ay (dXi / dt) + ωy2 .Xi D:1 -»n] wy = facteur de pondération Q = facteur de correction, du fait des marges d'erreurs possibles dans les données expérimentales, non indispensable.
La définition des autres paramètres et variables est la même que précédemment. Xi ; dXj / dt ; d2 Xj / dt2 et Xj sont des données fournies expérimentalement ou directement calculés à partir de ces données. Par exemple, ces données peuvent être obtenues par des criblages d'expressions d'ARNm.
Les inconnues de cette équation en sont les paramètres ( mi ; λy ; ω ; wy ), qui sont à fixer pour entièrement définir la relation en vue de simulations. On pose que quel que soit i et quel que soit j, | λy | > 1. .
Dans un autre aspect de l'invention, la relation entre Xj et l'ensemble des Xj est définie par une somme dont la pondération inclut un terme variable au cours du temps et tenant compte explicitement des vitesses respectives des modifications des Xj, vitesses représentées par les dérivées : dXj(t)/dt, noté dans la suite dXj/dt. Ceci revient à considérer que les vitesses de variation des Xj (des molécules j) influencent la résultante globale de leurs actions sur la cinétique de Xi (de la molécule i). On définit : (4) ay = ( dXj / dt ) / ∑ ( dXj / dt ) ; ay étant un facteur de pondération. D:1 ->n] ay est directement calculable, à partir des données quantitatives expérimentales, pour chaque temps expérimental. Ce facteur de pondération varie en fonction du temps. Dans ce cas, la relation entre Xi et l'ensemble des Xj est définie comme suit : (5) ∑ (ay . wy . Xj) + Q = mj .(d2Xi / dt2) + 2 .λy .(dXi / dt) + ωy 2 .Xi D:1->n]
Ou
(6) ∑ (ay . wy . Xj) + ci / df ) + 2 .λy (dXi / dt) + ωy .Xj D-1->n] La définition des paramètres et variables est la même que précédemment. Les équations (2) et (3) reviennent à un cas particulier des équations (5) et (6) où l'on pose arbitrairement : an = 1.
Dans un autre aspect de l'invention, la relation entre Xj et l'ensemble des Xj est définie par une somme dont la pondération inclut un terme variable au cours du temps et tenant compte explicitement des accélérations respectives des modifications des Xj, accélérations représentées par les dérivées secondes : d2Xj(t)/dt2 notées dans la suite d2Xj(t)/dt2. Ceci revient à considérer que les accélérations des variations des Xj (des molécules j) influencent la résultante globale de leurs actions sur la cinétique de Xi (de la molécule i). On définit : (7) ay = ( d2Xj / dt2 ) / ∑ ( û / dt2 ) ; ay étant un facteur de pondération.
Ici aussi, ay est directement calculable, à partir des données quantitatives expérimentales, pour chaque temps expérimental. Ce facteur de pondération varie en fonction du temps. Dans ce cas, la relation entre Xj et l'ensemble des Xj est définie comme précédemment : (8) ∑ (a . wy . Xj) + q = mi .(d2Xi / dt2) + 2 .λy .(dXi / dt) + ωy2.Xi D:1-»n]
Ou (9) ∑ (ay . Wy . Xj) + q =(d2Xj / dt2) + 2 .λy (dX| / dt) + ωy 2 .Xi [i:1-»n]
La définition des paramètres et variables est la même que précédemment. Seule la définition des paramètres ay (et par conséquent leur valeur et les valeurs calculées des autres paramètres) diffère par rapport à l'aspect précédent de l'invention. Ici aussi, Les équations (2) et (3) reviennent à un cas particulier des équations (8) et (9) où l'on pose arbitrairement : ay = 1.
Dans la mesure où les temps auxquels sont effectuées les mesures expérimentales sont longs par rapport aux accélérations, il est préférable de définir les pondérations par rapport aux vitesses que par rapport aux accélérations. - (2) Dans un autre aspect de l'invention, la relation entre Xi et les Xj est établie par une relation sigmoïde comportant un facteur de retardement associée à une relation de décroissance linéaire.
Dans ce cas, dans une mise en œuvre préférentielle, la relation entre Xj et
Figure imgf000066_0001
(10) (dX dt) = Kn . [ 1 / (1 + e'∑ wij Xj " bi) ] - K2i . Xi ; pour l'ensemble des sommets j ayant une action sur i (notés sommets j).
Dans cette formulation, la relation associée aux arrêtes correspondant à une interaction moléculaire inhibitrice pourra être l'inverse de celle associée aux arrêtes correspondant à une interaction activatrice (courbe sigmoïde décroissante ou croissante, respectivement), cette caractéristique étant directement obtenue lors du calcul des paramètres, par leurs signes positifs ou négatifs.
Le terme sigmoïde : Kιι . [ 1 / (1 + e-∑ w0J -bl) ] correspond au terme inertiel
Le terme :
Figure imgf000066_0002
correspond au terme de retour à l'état initial.
Xi = variable associée au sommet i. Xj = variable associée au sommet j. wij = facteur de pondération de l'effet de j sur i lorsque plusieurs sommets (j-i. J2, • Jn) agissent sur le sommet i. bi = facteur de retardement
K-π = facteur de limite maximale de variation de Xj . Son utilisation est liée au fait que le terme sigmoïde varie entre 0 (variation de 0 %) et 1 (variation de 100%). Ce facteur est donc requis pour rendre compte des situations ou Xi varie de plus de 100%. K2j = facteur de retour à l'équilibre.
Dans cet aspect de l'invention, la résultante pondérée de la combinaison des effets de l'ensemble des Xj sur Xj est incluse d'emblée. L'effet combiné de l'ensemble des Xj est représenté ici aussi par une somme pondérée, dans le terme inertiel. Dans un autre aspect de l'invention, la relation entre les Xj et les Xj est établie si ilairement à l'aspect précédent, mais avec l'introduction d'un facteur de pondération supplémentaire ay tel que défini par les équations (4) ou (7) précédentes. Dans ce cas, dans une mise en œuvre préférentielle, la relation entre X, et les Xj sera : (11) (dX dt) = Kn . [1 / (1 + e-∑aij' wij Xj " bi ) ] - K2i . Xi ; pour l'ensemble des sommets j ayant une action sur i (notés sommets j).
La définition des paramètres et variables est la même que précédemment ; le paramètre ay est défini soit par l'équation (4), soit par l'équation (7) précédentes. L'équation (10) revient à un cas particulier de l'équation (11) où l'on pose arbitrairement : a = 1.
Le paramètre ay est un terme variable au cours du temps et tenant compte explicitement des vitesses respectives des modifications des Xj, ou des accélérations respectives des Xj, selon que a est défini par l'équation (4) ou l'équation (7), respectivement. Dans la mesure où les temps auxquels sont effectuées les mesures expérimentales sont longs par rapport aux accélérations, il est préférable de définir les facteurs de pondération ay par rapport aux vitesses (équation
(4)) que par rapport aux accélérations (équation (7)).
Equation (4) :
(4) ay = ( dXj/ dt ) / ∑ ( dXj/ dt ) D-1-»n]
Equation (7) :
(7) ay = ( d2Xj / dt2 ) / ∑ ( d2 Xj/ dt2 ) D:1-»n] - (3) Cependant, d'autres types de relations mathématiques non citées ici comportant les caractéristiques décrites dans l'invention pourraient aussi être utilisées, leur utilisation dans le cadre de l'invention étant alors considérée comme tombant sous le coup du présent brevet. - (4) Il doit être entendu ici que l'utilisation de relations mathématiques non explicitées ici mais qui permettent de reprendre tout ou partie des caractéristiques décrites ci-dessus entrent dans le cadre de l'invention. Ainsi, il est possible, de façon non limitative, d'établir une relation entre Xi et Xj qui respecte l'existence d'inflexion de la courbe de Xi en fonction du temps, et une limite maximale de Xj dans l'intervalle des données expérimentales utilisées et dans les intervalles de temps expérimentaux par une fonction polynôme, ou une fonction sinus ou cosinus, etc.
Exemple 2 : Mise en oeuyre pratique de l'étape B)
De multiples techniques de procédures d'apprentissage, dont par descente de gradient au sein de graphes, sont disponibles pour effectuer le calcul des paramètres dans le domaine public (notamment en utilisant les transformées de Laplace, ou encore la méthode développée par Pearlmutter, 'Gradient calculations for dynamic récurrent neural networks : a suiNey', IEEE transactions on neural networks, 1995). Un autre exemple de méthode de calcul consiste à mette en œuvre la méthode de résolution numérique adaptative d'ordre S de Runge-Kutta (permettant d'utiliser un pas de temps non constant dans les données expérimentales) associée à un apprentissage par BPTT (back propagation through time). Le choix de la procédure d'apprentissage est notamment lié aux algorithmes utilisés pour définir la relation entre les couples (Xj, Xj). La personne de l'art pourra facilement effectuer ce choix, et le mettre en œuvre.
Le choix et la mise en œuvre d'une fonction d'erreur ne posent pas de difficulté particulière. En effet, plusieurs fonctions d'erreur peuvent être utilisées, et sont disponibles dans la littérature. A titre d'exemples, des types de fonctions d'erreur utilisables sont : E = ∑I[Xii(t)-X2i(t)]2.dt i t avec : Xij(t) : valeur de Xj au temps t calculée par simulation avec : X2j(t) : valeur de X; au temps t mesurée expérimentalement.
Ou encore l'erreur globale relative aux trajectoires données pour l'apprentissage : [∑∑(Xii(t)-X2i(t))2/∑∑(X2i(t))2] i t i t Λ/ = racine carrée du terme entre crochets [ ]. Xii(t) et X2ι(t ) étant définis comme ci-dessus. On peut également, pour calculer l'erreur relative locale (au niveau de la trajectoire d'un sommet du graphe), utiliser la formule : V[∑(X1i(t)-X2j(t))2/∑(X2i(t))2]
V = racine carrée du terme entre crochets : [ ]. Xϋ(t) et X2i(t ) étant définis comme ci-dessus.
Dans le cas de la mise en œuvre d'une relation Xj, Xj de type inertielle adaptée de l'oscillateur harmonique, et dans une mise en œuvre préférentielle, les contraintes suivantes seront imposées lors du calcul des paramètres :
Un seuillage sera introduit en imposant que pour toute relation élémentaire (Xi, Xj), les paramètres calculés respectent : λy > ωy (ce qui revient à imposer un amortissement important), ou encore, mi = valeur maximale calculée pour l'ensemble des relations (Xi, xj), ces deux critères pouvant être associés.
Au terme du calcul des paramètres, le graphe (ou réseau) est entièrement déterminé par les relations mathématiques associées aux arrêtes du graphe : il s'agit d'un réseau entièrement déterministe. Le graphe correspondant est orienté. Le réseau est peut être représenté de façon explicite par la mise en œuvre de techniques de représentation de réseaux de neurones utilisées en intelligence artificielle. Il s'agit d'un réseau non booléen, ni bayésien, ni organisé en couches, permettant de représenter des redondances des circuits et des boucles de rétro-action. Ce réseau déterministe permet la mise en œuvre de simulations sans coût de calcul notable, même pour un graphe de très grande taille.
Exemple 3 ; Mise en œuyre pratique de l'étape D)
Propagation :
De nombreuses méthodes de propagation sont disponibles dans la littérature, adaptées des technologies de réseaux de neurones développées en intelligence artificielle, et leur mise en œuvre ne pose pas de difficulté particulière à l'homme de l'art, le graphe dynamique étant entièrement déterministe à ce stade de la mise en œuvre.
A titre d'exemples de méthodes de propagation, on peut citer le logiciel Neural Network Toolbox 4.0.2, développé dans l'environnement de calcul
MATLAB, disponible à l'adresse internet http://www.mathtools.net/MATLAB/Neural Networks/index. html, commercialisé par la société The MathWorks, Inc. Ce logiciel de mise en oeuvre de réseau de neurones permet notamment de réaliser des propagations dans un réseau. D'autres exemples de propagations sont intégrés aux méthodes de Runge Kutta et de Pearlmutter citées dans le présent texte.
La propagation est inhérente aux méthodes d'apprentissage citées à l'exemple 2. L'étape de simulation consiste donc à utiliser la méthode de propagation mise en œuvre de l'étape B, ou toute autre méthode de propagation jugée adéquate par l'homme de l'art.
Plusieurs principes sont toutefois préférentiellement respectés lors de la mise en œuvre des simulations :
Dans une mise en œuvre particulière de l'invention, des procédures de seuillage sont associées à la méthode de propagation choisie, afin de diminuer les divergences (donc d'améliorer les convergences, c'est à dire la fiabilité). Celles-ci peuvent porter sur :
- Un seuillage inférieur (c'est à dire que toute valeur d'une variable prédite en dessous de 0 est ramenée à 0 (où à une valeur de bruit de fond minimale si celui-ci peut être défini par les données expérimentales) :
Seuillage : pour toute molécule i, quelle que soit la valeur Xi , Xj < 0 => Xj = 0.
- Un seuillage supérieur : en imposant un seuil maximal aux valeurs des Xj pouvant être obtenues lors des simulations (par exemple en fixant un seuil maximal correspondant à un facteur multiplicatif (pouvant de façon non exclusive être fixé entre 1 et 10) des valeurs maximales observées expérimentalement pour chaque Xj ; ce facteur peut éventuellement être défini lors de simulations de résultats expérimentaux réels disponibles en testant plusieurs valeurs de ce facteur. - L'introduction de contraintes dans les boucles lors des simulations : Ceci peut être réalisé par plusieurs méthodes, non exclusives les unes des autres, cette liste n'étant pas limitative. Toutes ces méthodes visent à imposer des contraintes, soit au nombre de boucles effectuées, soit aux gammes de valeurs des Xj : o Limitation du nombre de boucles pouvant être effectuées lors des simulations, le nombre maximal de boucles pouvant être défini à partir de l'analyse de données expérimentales réelles en tenant compte de la durée des simulations. o Utilisation du seuillage tel que décrit ci dessus pour éviter des phénomènes d' « explosion » dans des boucles, celles- ci allant alors se stabiliser au niveau maximal autorisé du seuil.
Itération :
La pratique d'itérations est couramment utilisée en informatique. Elle consiste ici à répéter la séquence de calculs de propagation en modifiant de façon systématique les stimuli. Elle peut ou non inclure une stratégie de calcul parallèle. Elle est facile à mettre en oeuvre par l'homme de l'art et ne nécessite pas d'être plus détaillée.
En permettant de décrire l'effet de telle ou telle modification du réseau sur l'ensemble du réseau, ces simulations visent à analyser le système biologique dans son ensemble, et donc à répondre aux enjeux cités plus haut. Ces simulations consistent donc à décrire l'évolution de l'ensemble des molécules du réseau d'interactions moléculaires au cours du temps suite aux stimuli « virtuels » initiaux, y compris, si cela s'avère biologiquement important, jusqu'à un nouvel état d'équilibre du graphe. Ces stimuli « virtuels » peuvent être appliqués de façon systématique sur chaque sommet du graphe, mais il est aussi possible d'effectuer plusieurs stimuli en même temps, ou de façon séquentielle, sur plusieurs sommets.
Proximité des états du graphe : Un calcul statistique de proximité entre chaque état final calculé et l'état voulu, ou entre chaque état final et l'état à modifier, est effectué. Il permet, pour chaque sommet, d'associer un critère statistique (proximité obtenue) au sommet sur lequel s'est exercé le stimulus et au stimulus exercé sur ce sommet.
De nombreuses possibilités de comparaison statistique des états d'un graphe existent, et leur choix ne pose pas de difficulté pour l'homme de l'art.
A titre d'exemple, la distance entre deux états d'un graphe peut être calculée comme suit :
Si Xii est la valeur de la variable Xj à l'état T du graphe ; Si X2i est la valeur de la variable Xi à l'état 2 du graphe, on peut calculer une distance mathématique entre les états 1 et 2 du graphe par : D = ∑ (Xii - X2i )2
Ou encore : D = / [ ∑ ( X1i - X2i )2 / ∑ ( X2j)2 ] i i Λ/ = racine carrée du terme entre crochets : [ ].
D'autres méthodes de statistiques classiques telles que comparaisons des populations, etc... sont aussi disponibles.
Les éléments de comparaison des graphes, sont de deux ordres :
- Points de convergence des valeurs Xj (end point), par exemple par comparaison des populations de Xj finaux entre états de graphes par statistiques classiques de comparaisons de populations (moyennes, variance...).
- Cinétiques de chaque molécule i (par exemple par comparaison des différences d'intégrales des courbes de cinétiques, et comparaison des populations de différences d'intégrales). Il est ainsi possible de comparer, par exemple, la cinétique des Xj au cours de l'établissement d'un processus pathologique à celle suite à un stimulus tel que défini plus haut, permettant de hiérarchiser les sommets et les stimuli par l'écart qu'ils provoquent par rapport au processus pathologique, en ne se limitant pas aux points de convergence. Dans ce cas, il est possible par exemple d'estimer la distance entre les deux cinétiques du graphe par les fonctions d'erreur citées plus haut : E = ∑ / [ Xii(t) - X2i (t)]2 . dt i t
Ou encore : V [ ∑ ∑ ( X1i(t) - X2i(t ) )2 / ∑ ∑ ( X2i(t) )2 ] i t i t
Ces deux types de comparaisons permettent d'évaluer statistiquement la proximité des graphes, dans les procédures de simulations décrites plus haut.
- D'autres types d'analyse d'analyse de proximité peuvent être mis en œuvre : voir par exemple : Pearlmutter, 'Gradient calculations for dynamic récurrent neural networks : a suπtey', IEEE transactions on neural networks, 1995. Le choix et la mise en œuvre de ces comparaisons sera facilement réalisée par l'homme de l'art et ne requiert pas plus de description ici.
Exemple 4 : modélisations et simulations à partir d'un graphe statique de 116 molécules Données biologiques utilisées :
1) Un graphe statique correspondant à un réseau d'interactions moléculaires dans la levure (Saccharomyces cerevisiae, organisme vivant eucaryote pouvant être considéré comme un système biologique répondant aux critères de complexité cités plus haut) a été construit par une saisie manuelle dans un fichier plat de type txt (sans mettre en œuvre de système automatisé particulier), à partir des données de la base de données KEGG : Kyoto Encyclopedia of Gènes and Génomes, données en accès libre à l'adresse internet http.7/www.αenome.ad.ip/keqq/keqq2.html.
Ce graphe est plus particulièrement centré sur les mécanismes de la respiration cellulaire (glycolyse, néoglucogénèse, métabolisme du Pyruvate et de l'acetyl CoA, etc.). Il comprend 116 molécules, enzymes ou facteurs de transcription. Il comprend 329 interactions uni- ou bi-directionnelles entre ces molécules.
Le graphe statique correspondant à ce réseau d'interactions moléculaires est donné ici à titre d'exemple, sous deux formes :
- Schéma (figure 2). Sur ce schéma, chaque rectangle représente .une protéine. Les lettres dans le rectangle sont les abréviations usuelles du nom de la protéine, selon la nomenclature de KEGG et de SGD : Saccharomyces Génome Database développée par le Department of Genetics at the School of Medicine, Université Stanford, USA.
- Tableau représentant les interactions moléculaires (Tableau 5 ci- dessous).
Ce tableau présente les données de graphe statique sous une forme directement utilisable par un système informatique de modélisation et de simulation tel que décrit dans l'invention.
Figure imgf000075_0001
Figure imgf000076_0001
Figure imgf000077_0001
Figure imgf000078_0001
Figure imgf000079_0001
Figure imgf000080_0001
Figure imgf000081_0001
Figure imgf000082_0001
Figure imgf000083_0001
Tableau 5 : Représentation du graphe sous forme de tableau
Le graphe représente les 329 interactions entre les 116 molécules du réseau. Les interactions sont représentées entre les molécules deux à deux.
Colonnes A : première molécule
Colonnes B : seconde molécule
Sens : sens de l'interaction : A-B : de A vers B B-A : de B vers A A<-> B : interaction dans les deux sens.
Ce tableau établit aussi la correspondance entre les codes ORF (open reading framë) de la base de données SGD (Dolinski, K., Balakrishnan, R., Christie, K. R., Costanzo, M. C, Dwight, S. S., Engel, S. R., Fisk, D. G., Hirschman, J. E., Hong, E. L., Issel-Taπ er, L., Sethuraman, A., Theesfeld,
C. L., Binkley, G., Lane, C, Schroeder, M., Dong, S., Weng, S., Andrada, R., Botstein, D., and Cherry, J. M. "Saccharomyces Génome Database " : http://www.yeastgenome.org/), et les abréviations des noms des protéines (elles aussi de SGD). Les codes ORF sont uniques pour une protéine donnée et permettent de l'identifier sans aucune ambiguïté. Ils permettent aussi d'établir un lien non ambigu avec les résultats de criblages sur puces à ADN (correspondance des séquences nucléiques des ARN messagers correspondants).
2) Des données de criblage d'expression d'ARN messagers sur puces à ADN concernant l'ensemble de ces gènes ont été saisies à partir de la publication : DeRisi JL, lyer VR, Brown PO : Exploring the metabolic and genetic control of gène expression on a genomic scale, Science. 1997 Oct 24;278(5338):680-6.
Cette publication décrit une expérience de culture de levures dans des conditions où la concentration de glucose dans le milieu de culture diminue progressivement (du fait de son utilisation par les levures pour la fermentation, du glucose n'étant rajouté à la culture à aucun temps de l'expérience). Au cours du temps, les levures présentent une modification de leur métabolisme, leur système respiratoire passant d'un fonctionnement en fermentation à un fonctionnement en respiration aérobie.
Cette culture de levures a été étudiée au cours du temps, notamment par la pratique de criblages d'expression de la quasi-totalité des ARN messagers de levure sur puces à ADN. Ces criblages ont été effectués à des temps successifs, les résultats produisant donc une cinétique de niveau d'expression pour chaque ARN messager. Les résultats montrent des variations du niveau d'expression d'un certain nombre d'ARN messagers au cours du temps, ceux-ci étant plus particulièrement nombreux parmi les ARN messagers des protéines de la respiration cellulaire, dont une partie importante est représentée dans le graphe décrit ci dessus. Dans ces conditions expérimentales, le graphe que nous avons construit présente donc un évolution dynamique au cours du temps, qui est représentée par les cinétiques des molécules du réseau (donc les sommets du graphe).
Il sera clair au lecteur que le graphe statique est déjà d'une taille trop grande, et comprend trop d'interactions et de boucles, pour permettre, même à un expert, de prédire correctement à partir du seul graphe statique son évolution dynamique telle qu'observée expérimentalement sans la mise en œuvre d'une méthode de modélisation dynamique adaptée.
L'ensemble des données expérimentales de criblage d'expression d'ARN messagers sur puces à ADN correspondant à cet article sont disponibles sur le site internet de l'Université de Stanford à l'adresse : http://cmgm.stanford.edu/pbrown/explore/array.txt.
A partir de ces données, les données expérimentales correspondant spécifiquement aux ARN messagers des molécules du graphe -ont été saisies manuellement (procédure de copier-coller dans un fichier plat de type txt) sous la forme suivante : chaque ligne correspond à une molécule du graphe, la première colonne identifiant l'ORF (open reading frame) par son code SGD, les colonnes suivantes correspondant aux mesures expérimentales. Les tableaux 6 à 8 ci-dessous donnent des exemples de données pour quelques molécules du graphe, données extraites de la page http://cmgm.stanford.edu/pbrown/explore/array.txt.
Figure imgf000085_0001
Tableau 6
Figure imgf000085_0002
Tableau 7
Figure imgf000085_0003
Tableau 8
Tableaux 6 à 8 : Exemples de données de criblages sur puces à ADN pour 2 des molécules du graphe, à partir des résultats de l'article : DeRisi JL, lyer VR, Brown PO : Exploring the metabolic and genetic control of gène expression on a genomic scale, Science. 1997 Oct 24;278(5338):680- 6. Les données complètes sont disponibles sur la page internet : http://cmgm.stanford.edu/pbrown/explore/array.txt. Compte tenu de la taille du tableau des données complètes, il n'en est montré ici qu'une partie. L'homme de l'art pourra très facilement récupérer les données correspondant aux autres molécules du graphe utilisé dans cet exemple, sur cette page internet qui correspond directement au tableau de l'ensemble des données.
Les criblages d'expression des ARN messagers ont été réalisés toutes les deux heures, pendant 12 heures, ce qui correspond à 7 temps expérimentaux (le temps initial plus les 6 temps suivants). Ceux-ci correspondent aux notations 1 à 7. Le lecteur trouvera toutes les explications correspondant à l'obtention de ces mesures dans l'article cité en référence.
Nom = abréviation du nom du gène (selon SGD)
ORF = code de l'open reading frame (selon SGD)
G = condition expérimentale correspondant à "l'état étalon" du graphe tel que décrit dans l'invention. Cet état étalon, dans cette série d'expériences, correspond à l'état initial de culture des levures. G1 , G2, G3, G4, G5, G6,
G7 correspondent tous au même échantillon biologique étalon.
R = états du graphe aux divers temps expérimentaux.
R1 correspond à l'état initial de culture des levures (même échantillon biologique que G1), au temps TO. R2 : TO + 2,5 heures, R3 : TO + 4 heures, R4 : TO + 6 heures, R5 : TO + 7,5 heures, R6 : TO + 9,5 heures, R7 : TO +
1 1 ,5 heures.
Les séries de valeurs G-Bkg et R-Bkg correspondent à des mesures absolues de signal. Par rapport au présent texte, les séries G-Bkg correspondent à xi0, et les séries R-Bkg correspondent à xit. G1 -Bkg = mesure de G1 moins le bruit de fond (background) lors des mesures expérimentales.
G2-Bkg = mesure de G2 moins le bruit de fond (background) lors des mesures expérimentales. G3-Bkg = mesure de G3 moins le bruit de fond (background) lors des mesures expérimentales.
G4-Bkg = mesure de G4 moins le bruit de fond (background) lors des mesures expérimentales.
G5-Bkg = mesure de G5 moins le bruit de fond (background) lors des mesures expérimentales.
G6-Bkg = mesure de G6 moins le bruit de fond (background) lors des mesures expérimentales.
G7-Bkg = mesure de G7 moins le bruit de fond (background) lors des mesures expérimentales. Les variations des valeurs mesurées sont liées aux variations des rendements des diverses réactions mises en oeuvre dans la méthode de mesure (puces à ADN) et justifient l'utilisation d'un état étalon comme référence de mesure.
R1-Bkg = mesure de R1 moins le bruit de fond (background) lors des mesures expérimentales.
R2-Bkg = mesure de R2 moins le bruit de fond (background) lors des mesures expérimentales.
R3-Bkg = mesure de R3 moins le bruit de fond (background) lors des mesures expérimentales.
R4-Bkg = mesure de R4 moins le bruit de fond (background) lors des mesures expérimentales.
R5-Bkg = mesure de R5 moins le bruit de fond (background) lors des mesures expérimentales. R6-Bkg = mesure de R6 moins le bruit de fond (background) lors des mesures expérimentales. R7-Bkg = mesure de R7 moins le bruit de fond (background) lors des mesures expérimentales.
R1. Ratio, R2. Ratio, R3.Ratio, R4.Ratio, Rδ.Ratio, Rβ.Ratio, R7.Ratio correspondent aux rapports R-Bkg / G-Bkg à chacun des 7 temps expérimentaux : ils correspondent aux variables telles que définies dans la description de l'invention : Xj = xit / xo.
On a donc produit deux tableaux sous la forme de fichiers plats de type txt, avec une correspondance mutuelle par le code ORF de chaque molécule.
Dans cet exemple de mise en œuvre, il n'a pas été nécessaire d'utiliser de système de base de données.
Mise en œuvre de la méthode :
Les étapes A et B ont été mises en œuvre comme suit :
Un lissage des données expérimentales a été effectué afin de disposer de plus de points temporels.
On a considéré que si la concentration en glucose dans le milieu de culture des levures avait été maintenue constante en supplémentant en glucose le milieu de culture, l'état initial, qui correspond aux mesures expérimentales au temps T0, serait un état stable. Ceci est en accord avec les données expérimentales disponibles et avec le texte de la publication dont ont été extraites les données.
Par ailleurs, on sait que si la culture de levures est à nouveau supplémentée en glucose après le temps 7 (T0 + 12 heures), elle va revenir à son état initial. Un état final du graphe a donc été défini comme suit : tendance du système biologique étudié à revenir à son état initial au temps To + 36 heures, et les données expérimentales obtenues à T0 ont été répliquées à ce temps. Ceci a été fait afin d'ajouter une contrainte, logique vis-à-vis de l'expérience, lors du calcul des paramètres des fonctions. Ceci n'est cependant pas à considérer comme une étape de mise en œuvre indispensable à la mise en œuvre de l'invention, mais comme un exemple de définition des états stables du système biologique à partir de cet exemple précis.
La relation utilisée entre les variables Xj correspondant à la molécule i et les variables Xj correspondant aux molécules j interagissant sur i a été la suivante :
(dXj/dt) = K-π . [1 / (1 + e-∑wij Xj-bi ) ] - K2i . Xi
Le calcul des paramètres a été effectué à partir des données expérimentales par une méthode classique d'apprentissage de réseaux de neurones, plus précisément à partir des algorithmes de la méthode de Runge Kutta de rétro-propagation dans le temps (BPTT : back propagation through time). Les calculs ont été effectués en double précision.
Le reste des méthodes mises en œuvre, qui ne posent pas de difficulté particulière à l'homme de l'art, ont été effectuées comme décrit plus haut.
Un graphe dynamique, entièrement déterministe a ainsi été obtenu, permettant de réaliser des simulations.
Résultats :
Lors des simulations, l'efficacité de la méthode a été vérifiée. En effet, le résultat de divergence moyen (erreur relative globale) des simulations par rapport aux données expérimentales est d'environ 0,30, cette divergence étant essentiellement concentrée sur 8 sommets (molécules) du graphe, pour cette série de données, alors que pour les 108 sommets (molécules) restant, la divergence est très faible. Ce résultat d'erreur relative globale montre que les cinétiques calculées lors des simulations sont proches des données réelles, car des cinétiques aléatoires auraient donné un calcul d'erreur supérieur à 1.
La divergence globale des simulations par rapport aux données expérimentales sur l'ensemble du graphe et l'ensemble des cinétiques a été estimée par le calcul d'erreur relative suivant :
Erreur globale relative: V [ ∑ ∑ ( Xii(t) - X2i(t ) )2 / ∑ ∑ ( X2i(t) )2 ] i t i t V = racine carrée du terme entre crochets [ ]. X-ii(t) = valeur de X calculée au temps t de la simulation, X2i(t ) = valeur de X mesurée expérimentalement au temps t. ∑ = somme des valeurs aux différents temps t
Ce résultat de simulation est satisfaisant, d'autant plus si l'on tient compte du taux d'erreurs de mesures, puisqu'il est légèrement inférieur au taux d'erreurs de mesures lors des expériences ayant servi à générer les données expérimentales sur puces à ADN.
Le taux de non-reproductibilité des données expérimentales peut être estimé par le rapport R1. Ratio des données expérimentales (Tableau 8), et est globalement de 14% dans cet exemple.
Ce résultat de divergence globale est obtenu par un calcul d'erreur relative permettant de comparer deux cinétiques (ou trajectoires) dans leur ensemble. Il ne peut naturellement pas être utilisé pour estimer le taux de non-reproductibilité des données expérimentales puisqu'on ne dispose ici que d'une seule cinétique expérimentale, pour ces conditions expérimentales, pour chaque molécule du réseau. Le calcul de ce taux de non-reproductibilité a donc été effectué par la moyenne des ratios R1. Le fait que ces deux calculs d'erreurs soient différents ne permet pas de les comparer directement au sens strict. Cependant, on voit que l'erreur relative globale des simulations et le taux de non-reproductibilité des mesures expérimentales sont proches : 0,3 et 0,14 respectivement, 1 étant le seuil au dessus duquel les simulations et les mesures peuvent être considérées comme non-fiables.
Bien qu'il soit possible par la méthode de l'invention de descendre lors des simulations à un résultat de divergence inférieur au taux de non- reproductibilité des données expérimentales, il est clairement inutile de descendre à un résultat de divergence inférieur "à cette limite de reproductibilité des données expérimentales, quelles que soient celles-ci. En effet, puisque des données expérimentales sont utilisées pour le calcul des paramètres, cela reviendrait à introduire tout de même un risque de divergence vis-à-vis du phénomène biologique réel étudié, dont la divergence vis-à-vis des mesures expérimentales peut être estimée égale au taux de non-reproductibilité des expériences de mesure, sans que l'on puisse prédire le sens de cette divergence (qui peut de plus varier en fonction des molécules du réseau).
A titre d'exemple, le tableau 9 donne le détail des calculs des divergences de l'ensemble des cinétiques lors des simulations pour l'ensemble des molécules du réseau, sous la forme d'un tableau récapitulatif.
Figure imgf000092_0001
pour l'ensemble des molécules du réseau Les divergences ont été estimées pour chaque molécule du réseau par le calcul de l'erreur relative sur l'ensemble de la trajectoire de la molécule concernée, suivant la formule suivante (erreur relative locale) : V [ ∑ ( X1i(t) - X2i(t ) )2 / ∑ ( X2i(t) )2 ]
Xii(t) = valeur de Xi calculée au temps t de la simulation, X2i(t ) = valeur de Xi mesurée expérimentalement au temps t, ∑ = somme des valeurs aux différents temps t = racine carrée du terme entre les crochets [ j -
Ce calcul revient à calculer la différence d'intégrale entre les courbes des cinétiques observées expérimentalement et les cinétiques calculées lors des simulations. Elle concerne donc aussi bien l'ensemble de la cinétique que l'état final.
Dans une variante de cet exemple, la modélisation et les simulations ont été mises en œuvre de la même manière que décrite ci-dessus, et à partir des mêmes données biologiques, avec la seule modification suivante lors du calcul des paramètres par rétro-propagation dans le temps :
Les variables associées aux sommets du graphe ne recevant pas d'arc ou arrête, c'est-à-dire correspondant aux molécules ne recevant pas d'interaction ("inputs" du graphe) ont été exclues du calcul d'erreur globale lors de l'apprentissage, leurs valeurs restant donc fixées aux valeurs expérimentales mesurées pendant cette procédure. Ceci a été effectué : - afin d'éviter de simuler des cinétiques sur ces sommets lors de la descente de gradient ce qui risque de majorer les erreurs, - et car ces sommets ne recevant pas eux mêmes d'inputs, leurs cinétiques sont de fait indépendantes des résultats de calculs des paramètres des relations mathématiques reliant les sommets. En d'autres termes, seules les molécules recevant au moins une interaction (arrête orientée vers le sommet du graphe leur correspondant) ont été prises en compte pour le calcul d'erreur lors de l'apprentissage.
Afin d'éviter toute confusion, cette variante ne consiste bien sûr pas à enlever du graphe les sommets "inputs", mais à imposer que leur cinétique reste la cinétique mesurée expérimentalement, ceci uniquement lors des simulations pratiquées pendant les calculs d'erreur de la procédure de calcul des paramètres par rétro-propagation dans le temps. Les paramètres des relations mathématiques reliant ces sommets à d'autres sommets du graphe sont donc calculés, comme pour toutes les autres arrêtes du graphe, et le modèle dynamique finalement obtenu inclut ces sommets.
Dans cette variante, les résultats de simulation obtenus ont été similaires à ceux montrés ci-dessus, bien que légèrement meilleurs.
La figure 3 donne à titre d'exemple les cinétiques mesurées expérimentalement et les cinétiques calculées par simulation pour quelques gènes représentatifs de l'ensemble des résultats obtenus par la mise en œuvre de cette variante.
Exemple 5 : Modélisations, simulations et validation de la capacité prédictive à partir d'un graphe statique de 133 molécules
Cet exemple montre la mise en œuvre de l'ensemble de la méthode (étapes A et B ou A', puis C, D, E) et son efficacité prédictive dans une application similaire à une identification / sélection de cibles thérapeutiques.
Données biologiques utilisées :
1) Un graphe statique correspondant à un réseau d'interactions moléculaires dans la levure (Saccharomyces cerevisiae, organisme vivant pouvant être considéré comme un système biologique répondant aux critères de complexité cités plus haut) a été construit selon les mêmes principes que dans l'exemple 4. Ce graphe inclut plus particulièrement le graphe de l'exemple 4, mais avec des molécules et des interactions supplémentaires. Il comprend 133 molécules, enzymes ou facteurs de transcription. Il comprend 407 interactions uni- ou bi-directionnelles entre ces molécules.
Le graphe statique correspondant à ce réseau d'interactions moléculaires est donné ici à titre d'exemple, sous deux formes : - Schéma (figure 4). Les principes de représentation et les commentaires explicatifs sont les mêmes que pour l'exemple 4 (figure 2). - Tableau représentant les interactions moléculaires additionnelles par rapport au graphe statique de l'exemple 4 (Tableau 10 ci-dessous). Le tableau complet représentant le graphe statique utilisé dans cet exemple 5 est donc l'addition des tableaux 5 et 10.
Figure imgf000095_0001
Figure imgf000096_0001
Figure imgf000097_0001
Tableau 10 : nteractions moléculaires additionnelles par rapport au graphe statique de l'exemple 4
Les commentaires du tableau 10 sont similaires à ceux du tableau 5.
2) Des données de criblage d'expression d'ARN messagers sur puces à ADN concernant l'ensemble de ces gènes ont été saisies à partir de la même publication que dans l'exemple 4 et selon les mêmes principes : DeRisi JL, lyer VR, Brown PO : Exploring the metabolic and genetic control of gène expression on a genomic scale, Science. 1997 Oct 24 ;278(5338) :680-6.
Par ailleurs, 3 métabolites ont été introduits parmi les 133 molécules du réseau. Contrairement aux autres molécules du réseau, ces métabolites n'ont naturellement pas d'ARN messager correspondant. Leurs valeurs ont été définies comme suit :
Glucose : ses concentrations au cours du temps ont été mesurées par les auteurs de la publication, aux mêmes temps que ceux auxquels ont été pratiqués les mesures d'expression d'ARN messagers. Les concentrations correspondantes sont données graphiquement dans la figure 4 de l'article cité. Afin d'exprimer ces valeurs sous forme de ratio, chaque valeur de la concentration en Glucose dans le milieu de culture à un temps donné a été divisée par la concentration en Glucose au temps initial de l'expérience (ceci afin de mesurer les ratios par rapport au même référentiel que pour les mesures d'ARN messagers, dont les ratios sont exprimés par le rapport de la mesure au temps t divisée par la mesure au temps initial). Cette variable associée au glucose correspond bien aux variables telles que définies dans la description de l'invention : Xi = Xjt / x0. en résulte pour le glucose les valeurs de variable suivantes
Figure imgf000098_0001
Tableau 11 : Valeurs de la variable associée à la molécule du graphe Glucose
Acétate et Glutamate : les concentrations de ces molécules n'ont pas été mesurées par les auteurs. Il a donc été décidé d'extrapoler des valeurs pour ces molécules à partir de la connaissance du système biologique étudié et de la description des expériences dans l'article. Dans la mesure où cette expérience est essentiellement fondée sur la chute progressive de la concentration en glucose dans le milieu de culture et où les autres paramètres du milieu de culture sont en première approximation considérés comme constants, il a été considéré que les concentrations du Glutamate et de l'Acétate, respectivement, étaient constantes au cours de l'expérience.
Le fait de travailler avec des ratios permet donc de fixer leurs valeurs selon les mêmes principes que pour le Glucose :
Xi = Xit / Xθ - Donc au temps initial (T0 = 0), Xj0 = xio / Xio = ,
Et, la valeur de Xj étant considérée comme constante au cours du temps, elle reste toujours égale à 1.
Il en découle le tableau de valeurs suivant
Figure imgf000099_0001
Tableau 12 : Valeurs des variables associées respectivement à la molécule du graphe : Glutamate et à la molécule du graphe : Acétate
On a donc produit deux tableaux sous la forme de fichiers plats de type txt, avec une correspondance mutuelle par le code ORF de chaque molécule (concernant les protéines / ARN messagers) ou le nom de molécule (concernant le Glucose, le Glutamate et l'Acétate). Dans cet exemple de mise en œuvre il n'a pas été nécessaire d'utiliser de base de données.
Mise en œuvre de la méthode :
Les étapes A et B ou l'étape A' ont été mises en œuvre de manière similaire à l'exemple 4 :
Un lissage des données a été effectué par extrapolation linéaire.
On a considéré que l'état initial, qui correspond aux mesures expérimentales au temps TO, serait un état stable si le milieu de culture avait été maintenu constant en le supplémentant en glucose. A la différence de l'exemple 4, on n'a pas défini un état final du graphe qui correspondrait à un retour à l'état initial suite à une supplémentation en glucose après le temps 7 (TO + 12 heures). Les seules données expérimentales à avoir été utilisées pour le calcul des paramètres ont donc été les données effectivement décrites dans l'article et présentes sur le site internet de l'Université de Stanford à l'adresse : http://cmqm.stanford.edu/pbrown/explore/array.txt, et correspondantes aux molécules du réseau, sans aucune extrapolation autre que celle concernant les molécules Glucose, Glutamate et Acétate décrites pjus haut.
De même que dans l'exemple 4, la relation utilisée entre les variables X, correspondant à la molécule i et les variables Xj correspondant aux molécules j interagissant sur i a été la suivante :
(dX dt) = Kϋ . [1 / (1 + e-∑wij Xj-bi ) ] - K2i . Xi
Le pas d'apprentissage pour le calcul des paramètres a été fixé à 16 minutes.
Le calcul des paramètres et le reste des méthodes mises en œuvre a été effectué comme décrit dans l'exemple 4, aboutissant à l'obtention d'un graphe dynamique, entièrement déterministe, permettant de réaliser des simulations.
L'étape C a été mise en œuvre comme suit :
L'objectif a été de montrer la capacité de la méthode à prédire un résultat nouveau non utilisé pour la construction du graphe dynamique. Toujours dans le même article : DeRisi JL, lyer VR, Brown PO : Exploring the metabolic and genetic control of gène expression on a genomic scale, Science. 1997 Oct 24 ;278(5338) :680-6, les auteurs ont aussi effectué un criblage d'expression d'ARN messagers d'une souche génétiquement modifiée de levure, dans laquelle a été effectué le « knock out » (« délétion ») du gène TUP1 (code de son ORF dans la base de données SGD : YCR084C), présent dans le graphe statique de 133 molécules. Ces données n'ont pas été utilisées pour la construction du graphe dynamique lors des étapes A et B.
Les conditions de culture et de criblage de cette souche sont amplement décrites dans l'article, mais pour plus de clarté dans la description de cet exemple de mise en œuvre, on peut noter les points suivants concernant ce criblage : la souche génétiquement modifiée a été cultivée dans les mêmes conditions de culture que la souche sauvage utilisée pour le reste des expériences, en présence de glucose, ce qui correspond pour la souche sauvage aux conditions de culture au temps initial des autres criblages effectués et décrits dans l'exemple 4 (T = T0). Ce criblage a été effectué en mesurant les rapports entre le niveau d'expression de chaque gène dans la souche présentant la délétion du gène TUP1 par rapport au niveau d'expression du même gène dans la souche ne présentant pas de délétion (souche sauvage). Ces données sont donc exprimées par rapport au même référentiel de mesure que celles décrivant les cinétiques lors de la privation de glucose (voir exemple 4), ce référentiel correspondant au temps initial des autres criblages effectués et décrits dans l'exemple 4 (T =
T0).
Afin de montrer la capacité de la méthode à sélectionner de façon pertinente des molécules-cibles sur lesquelles une action biologique ou pharmacologique permet de faire évoluer le système biologique étudié vers un état donné, on a donc utilisé le graphe dynamique obtenu par la mise en œuvre des étapes A et B pour poser la question suivante : « où faudrait-il agir sur le réseau de 133 molécules pour faire évoluer ce réseau vers un état le plus proche possible de l'état décrit par le criblage d'expression de la souche présentant la délétion du gène TUP1 ? » Cette question est 5 exactement du même type que celles posées dans la description de la mise en œuvre de la méthode en vue de la sélection de cibles thérapeutiques.
Dans la mesure où la souche de levure présentant la délétion du gène TUP1 ne diffère initialement de la souche « sauvage » que par cette
10 délétion, on a donc défini « l'état à modifier » du graphe comme étant l'état de référence de la souche sauvage cultivée en présence de glucose, c'est à dire son état au temps initial des autres criblages effectués, décrits dans l'exemple 4 (T = T0) et dont les résultats pour 130 des molécules du réseau (autres que Glucose, Glutamate et Acétate) sont disponibles sur le site
15 internet de l'Université de Stanford à l'adresse : http://cmgm.stanford.edu/pbrown/explore/array.txt.
Les données de criblage d'expression d'ARN messagers correspondant à cet état ont donc été celles du temps initial utilisé pour la construction du
20 graphe dynamique. Ces données permettent de définir numériquement l'état à modifier : une valeur numérique de ratio d'expression est associée à chaque molécule du réseau ; concernant les trois métabolites Glucose, Glutamate et Acétate, leur valeur dans l'état à modifier a aussi été leur valeur au temps 0 telle que décrite plus haut. Cette définition de l'état à
25 modifier est donnée ici à titre d'exemple. Un autre état à modifier aurait pu être défini par l'homme de l'art face à la mise en œuvre de l'invention pour d'autres applications.
L'étape D a été mise en œuvre comme suit
J o Cette étape consiste à pratiquer des simulations itératives, telles que décrites dans l'invention.
La question posée à laquelle les simulations devaient répondre a été : « où faudrait-il agir sur le réseau de 133 molécules pour faire évoluer ce réseau vers un état le plus proche possible de l'état décrit par le criblage d'expression de la souche présentant la délétion du gène TUP1 ? »
La souche de levure présentant la délétion du gène TUP1 ne diffère initialement de la souche « sauvage » que par cette délétion. Cette délétion revenant à une inhibition constante de l'expression du gène TUP1 , les simulations ont consisté à simuler, de façon itérative, l'inhibition constante de chacune des 133 molécules du réseau, et à effectuer un calcul de propagation au cours du temps de cette inhibition au sein du réseau. Pour chaque simulation, une seule molécule du réseau a été inhibée, puisque l'état à atteindre correspond à une évolution du système biologique modélisé (la souche de levure) suite à une seule inhibition (délétion du gène TUP1). On a donc réalisé 133 simulations.
D'après les commentaires des auteurs, les données expérimentales de l'article et les données de criblage d'expression d'ARN messagers concernant la souche présentant la délétion du gène TUP1 (accessibles sur le site internet de l'Université de Stanford à l'adresse : http://cmgm.stanford.edu/pbrown/explore/tupsearch.html). la délétion du gène était incomplète dans cette expérience biologique, le ratio : [niveau d'expression du gène TUP1 dans la souche délétée] / [niveau d'expression du gène TUP1 dans la souche sauvage] étant égal à 0,1 dans une mesure, et à 0,45 lors de la réplication de la mesure (moyenne : 0,28) . Dans le cas d'une délétion complète ce ratio aurait été égal en théorie à 0, et égal en pratique au bruit expérimental de mesure. Pour la mise en œuvre des simulations, afin de pouvoir reproduire une inhibition de type délétion, on a donc défini numériquement l'inhibition, pour chaque molécule du graphe, comme la multiplication par un facteur 0,1 du niveau d'expression de cette molécule au temps initial (état à modifier tel que défini plus haut), ce facteur correspondant à la valeur de l'inhibition la plus forte mesurée expérimentalement pour ce gène.
Donc pour chacune des 133 simulations effectuées, la simulation a consisté à imposer une valeur X-, constante dans le temps telle que X-, = [0,1 . XJO ] à une molécule unique i du graphe (X,o = valeur, d'expression (ratio) de la molécule i au temps expérimental initial T = T0), les valeurs des Xi des autres molécules étant initialement fixées à leur valeur dans l'état à modifier défini plus haut, et libres d'évoluer dans le graphe dynamique en fonction des calculs de propagation. Pour chacune des 133 simulations, la molécule i a été différente : l'effet de chaque inhibition d'un facteur 0,1 sur chaque molécule du graphe a été testé de façon systématique.
Dans cet exemple de mise en œuvre, l'inhibition a été imposée comme constante dans le temps lors des simulations : ainsi, lors des calculs de propagation, un éventuel retour de propagation sur la molécule i inhibée
(« feedback ») a été sans effet sur cette inhibition (Xj restant stable à sa valeur initiale de simulation). Ceci a été effectué afin de reproduire l'effet de la délétion d'un gène, qui est elle-même constante dans le temps. Ceci n'est cependant pas un pré-requis de la mise en œuvre de l'invention. Dans la mise en oeuvre de l'invention pour d'autres applications, l'homme de l'art peut décider de simuler des activations ou des inhibitions non constantes dans le temps, où à des temps différents.
Le calcul de propagation au sein du graphe suite à chacune des 133 inhibitions a été poursuivi pendant une durée simulée de 12 heures. Ces éléments étant posés, l'étape D a été mise en œuvre comme décrit dans l'invention, sans particularité notable, et sans présenter de difficulté particulière pour l'homme de l'art. Les calculs de simulations, consistant à propager l'inhibition initiale au cours du temps ont été réalisées par les mêmes principes et les mêmes outils que les simulations faisant partie de la procédure de calcul des paramètres.
Chacune des 133 simulations de l'étape D a ainsi aboutit, au temps 12 heures (durée de la propagation simulée), au calcul d'une nouvelle valeur numérique associée à chaque molécule du réseau, définissant un état du graphe : « état obtenu par simulation ». On a donc obtenu 133 « états obtenus par simulation » différents.
L'étape E a été mise en œuvre comme suit :
Cette étape consiste à hiérarchiser les molécules du graphe, et les effets exercés sur ces molécules lors des simulations, en référence à la proximité plus ou moins grande de la résultante de ces effets avec un état du graphe à atteindre.
Dans cet exemple de mise en œuvre, l'état à atteindre a été l'état de la souche de levure présentant une délétion du gène TUP1 décrite plus haut.
Ces données de criblage d'expression d'ARN messagers dans les conditions de délétion du gène TUP1 sont disponibles sur le site internet de l'Université de Stanford à l'adresse : http://cmgm.stanford.edu/pbrown/explore/tupsearch.html . Elles ont été saisies manuellement pour chaque molécule du graphe par une requête concernant l'ORF de cette molécule à cette adresse et insérées dans un fichier plat de type txt sous la forme (par exemple pour le gène CIT1 ) :
Figure imgf000106_0001
Tableau 13 : Exemple de donnée de criblage sur puce à ADN pour une molécule du graphe, dans la condition expérimentale de délétion du gène TUP1, à partir des résultats de l'article : DeRisi JL, lyer VR,- Brown PO : Exploring the metabolic and genetic control of gène expression on a genomic scale, Science. 1997 Oct 24 ;278(5338) :680-6. Pour chaque molécule i, la valeur du Xj mesurée expérimentalement correspond à la colonne « Avg. R/G » dans les données expérimentales accessibles sur le site internet : http://cmgm.stanford.edu/pbrown/explore/tupsearch.html.
La valeur du Xj correspondant au Glucose pour la souche présentant la délétion du gène TUP1 a été fixée à 1 puisque le criblage a été pratiqué sur une culture en présence de glucose. Les valeurs des Xj correspondant au Glutamate et à l'Acétate pour la souche présentant la délétion du gène TUP1 ont été fixées à 1 puisque le criblage a été pratiqué sur une souche dans un milieu de culture identique à celui de la souche sauvage au temps 0 (entre les deux cultures, les ratios des métabolites dans le milieu de culture sont donc égaux à 1).
Figure imgf000106_0002
Figure imgf000107_0001
Tableau 14 : Liste complète des valeurs de l'état à atteindre tel que défini plus haut
L'ensemble de ces valeurs définit donc numériquement un état du graphe « l'état à atteindre ». L'étape E consiste alors à calculer la distance entre d'une part « l'état à atteindre » du graphe, et d'autre part chacun des 133 « états obtenus par simulation » du graphe obtenus à l'étape D.
Ce calcul de distance est décrit précédemment (proximité des états du graphe) et ne pose pas de difficulté particulière à l'homme de l'art. Il consiste à comparer deux états du graphe en comparant deux à deux l'ensemble des valeurs Xi associées à chaque molécule i du graphe.
Dans cet exemple précis, le calcul de distance utilisé a été effectué en deux étapes :
L'étape 1 a consisté à effectuer une première classification par des calculs de distance de façon classique :
Distance d'ordre 1 : somme des valeurs absolues des différences entre les valeurs des Xj mesurées expérimentalement lors de la délétion du gène TUP1 (X,2 dans la formule ci-dessous) et les valeurs des Xj mesurées par simulation (Xπ dans la formule ci-dessous)
Figure imgf000108_0001
i
On a donc obtenu 133 calculs de distance, chacun correspondant à la distance entre d'une part l'état obtenu par simulation d'une propagation de 12 heures suite à l'inhibition d'une des molécules du graphe d'un facteur
0,1 et d'autre part l'état à atteindre.
Ces 133 distances calculées ont ensuite été classées en ordre de valeur croissant (de la plus grande proximité avec l'état à atteindre vers la plus grande distance avec l'état à atteindre). Cette classification correspond directement à la classification des molécules du graphe, de celle dont l'inhibition fait évoluer le graphe vers un état le plus proche de l'état à atteindre, à celle dont l'inhibition fait évoluer le graphe vers un état le plus éloigné de l'état à atteindre : il en a résulté une classification directe, et donc une hiérarchisation, des molécules-cibles sur lesquelles agir par inhibition pour faire évoluer le graphe vers l'état qu'il présente lorsque le gène TUP1 est délété.
L'étape 2 a consisté, à la suite de cette première classification, à l'affiner. Les molécules mieux classées que les molécules « outputs » du graphe lors du classement précédent (distances les plus faibles) ont été classées à nouveau entre elles par un second calcul de distance plus qualitatif : la différence entre la « sensibilité » et la « spécificité » des simulations : distance = sensibilité - spécificité. Cette étape de classification est donnée ici à titre d'exemple mais n'est pas indispensable à la mise en œuvre de l'étape E.
La sensibilité et la spécificité des simulations ont été calculées comme suit :
A partir des données expérimentales mesurées lors du criblage d'expression d'ARN messagers de la souche de levure présentant une délétion du gène TUP1 , on a identifié toutes les molécules du graphe présentant une variation d'expression supérieure à un facteur 2 par rapport à l'état de référence (souche de levure sauvage au temps initial T = T0 en présence de glucose), soit un groupe A de molécules.
De même, pour chaque « état obtenu par simulation » on a identifié toutes les molécules du graphe présentant une variation d'expression supérieure à un facteur 2 par rapport à l'état de référence (souche de levure sauvage au temps initial T = T0 en présence de glucose), soit un groupe Bj de molécules. Bi = groupe de toutes les molécules du graphe présentant une variation d'expression supérieure à un facteur 2 par rapport à l'état de référence suite à la simulation de l'inhibition de la molécule i du graphe.
La sensibilité a alors été définie, pour chacune des 133 simulations, comme le nombre de molécules du groupe Bj effectivement présentes dans le groupe A. Cela revient à évaluer, pour les variations quantitativement importantes d'expression des molécules (supérieures à un facteur 2) si la simulation induit effectivement les variations présentes dans les données expérimentales de l'état à atteindre. Plus la valeur de la sensibilité est élevée, plus la distance entre les deux états du graphe comparés est faible.
La spécificité a alors été définie, pour chacune des 133 simulations, comme le nombre de molécules du groupe Bi absentes du groupe A. Cela revient à évaluer, pour les variations quantitativement importantes d'expression des molécules (supérieures à un facteur 2) si la simulation n'induit pas des variations d'expression absentes dans les données expérimentales de l'état à atteindre. Plus la valeur de la spécificité est faible, plus la distance entre les deux états du graphe comparés est faible.
La différence sensibilité - spécificité revient donc à évaluer la distance sur les critères combinés de l'induction par la simulation des variations d'expression présentes dans l'état à atteindre et de la non-induction par la simulation de variations d'expression absentes dans l'état à atteindre.
Ces deux calculs (sensibilité et spécificité) reviennent simplement à compter pour chaque état du graphe le nombre de variables Xi dont la valeur est supérieure à 2 ou inférieure à 0,5 et ne posent aucune difficulté à l'homme de l'art. Ils peuvent d'ailleurs être effectués manuellement.
La différence entre les deux valeurs, sensibilité - spécificité, est elle aussi très simple et peut par exemple être calculée de façon manuelle, ou automatique par un tableau de logiciel Excel (Microsoft) ou tout autre logiciel équivalent.
Résultats :
Résultats de la mise en œuvre des étapes A et B ou de l'étape A' :
Lors des simulations, l'efficacité de la méthode a été vérifiée.
Dans l'exemple montré ici, le calcul d'erreur- globale relative d'apprentissage a été de 25,90 %, ce qui est satisfaisant. Ce résultat d'erreur relative globale montre que les cinétiques calculées lors des simulations sont proches des données réelles ; des cinétiques aléatoires auraient donné un calcul d'erreur supérieur à 1.
Exemple de courbes de paramétrage
La figure 5 donne à titre d'exemple les cinétiques mesurées expérimentalement (en blanc) et les cinétiques calculées par simulation (en noir) pour quelques molécules représentatives de l'ensemble des résultats obtenus par la mise en œuvre de cette variante des étapes A et B ou de l'étape A'.
On voit que ce résultat est très satisfaisant, d'autant plus si l'on tient compte des erreurs de mesures. Les considérations de l'exemple 4 à ce sujet restent pertinentes ici aussi. Le calcul des paramètres effectué à l'étape B a donc permis d'obtenir un graphe dynamique rendant bien compte des données expérimentales utilisées pour le calcul.
Résultats de la mise en œuvre des étapes C, D et E (capacité prédictive de l'invention) : Ces trois étapes aboutissent à la classification hiérarchique des molécules par classification hiérarchique des distances calculées entre d'une part l'état à atteindre (délétion du gène TUP1) et d'autre part les 133 états obtenus par simulation.
Lors de la mise en œuvre de ces étapes de simulations, l'efficacité de la méthode a été vérifiée.
Le résultat de la première étape de classification est résumé à titre d'exemple dans le tableau suivant, sous la forme d'un tableau récapitulatif. Chaque molécule du réseau est classée par la distance entre l'état à atteindre du graphe et l'état du graphe obtenu par la simulation de l'inhibition constante par un facteur 0,1 de cette molécule.
Figure imgf000112_0001
Figure imgf000113_0001
Tableau 16 : Distances entre l'état à atteindre du graphe et l'état du graphe obtenu par la simulation de l'inhibition constante par un facteur 0,1 de chaque molécule
Le classement par ordre croissant des distances donne directement le classement des molécules de celle dont l'inhibition est la plus susceptible de provoquer l'état à atteindre du graphe à celle dont l'inhibition est la moins susceptible de la provoquer. On voit que la molécule TUP1 est classée en première position, ce qui est bien le résultat attendu. La méthode est donc validée.
La figure 6 donne à titre d'exemple une représentation schématique de ce résultat de classification des molécules du réseau. Chaque point correspond à une molécule du réseau. Les ordonnées correspondent aux valeurs de distance calculées. En abscisse les 133 molécules du réseau sont classées de gauche à droite de celle associée à la distance la plus faible à celle associée à la distance la plus élevée.
Il est évident qu'on a bien obtenu directement une classification des molécules du graphe. Celle-ci est de même nature et a été obtenue selon les mêmes méthodes que celle qui serait obtenue dans une application de l'invention pour la recherche de cibles thérapeutiques.
Selon cette classification, la molécule TUP1 est classée en première position. Ce résultat est tout à fait satisfaisant. A titre d'exemple, le test expérimental des 5 premières molécules telles que classées ici donnerait donc un taux de succès de 100% pour la sélection de la molécule pertinente.
Cette classification présente aussi l'intérêt de pouvoir définir un « bornage » de l'ensemble des molécules-cibles sélectionnées. En effet, certaines molécules du graphe n'envoient pas d'interaction vers une autre molécule (ce sont donc des « sorties » ou « outputs » du graphe : « molécules outputs » dans la suite du texte) ; la simulation de l'inhibition de ces molécules n'entraîne donc pas de propagation de l'inhibition au sein du graphe qui reste donc globalement stable (puisqu'on a considéré que l'état initial était stable). De ce fait, ces molécules sont d'ailleurs classées dans un groupe contigu, de la 27ième position à la 30ième position. I_es molécules moins bien classées que les molécules outputs ne sont donc pas intéressantes en tant que molécule-cible.
En d'autres termes, le classement peut être interprété comme suit : la molécule la mieux placée est celle dont la simulation d'inhibition aboutit à l'état du graphe le plus proche de l'état à atteindre. Pour les molécules suivantes on s'éloigne progressivement de l'état à atteindre, en se dirigeant vers l'état à modifier qu'on atteint lorsqu'on arrive aux molécules outputs (cet état n'étant pas modifié lors de la simulation de l'inhibition des molécules outputs, à la molécule output près). Au delà des molécules outputs, les simulations d'inhibition aboutissent à des états du graphe qui s'éloignent progressivement à la fois de l'état à atteindre et de l'état à modifier. On a bien aboutit à la sélection d'un nombre limité de molécules-cibles (celles qui sont mieux classées que les outputs, ici 26 molécules), qui sont elles-mêmes hiérarchisées en terme de priorité. Le classement de la molécule TUP1 montre que cette hiérarchisation est satisfaisante. L'étape 2 de classification (le calcul sensibilité - spécificité), bien que non indispensable compte-tenu du résultat qui précède, a été ensuite mise en œuvre afin d'améliorer la classification des cibles. Elle n'a été appliquée qu'aux molécules du graphe mieux classées que les molécules outputs à l'étape précédente (26 molécules). Le classement ainsi obtenu est donné dans le tableau 17.
Figure imgf000115_0001
Figure imgf000116_0001
Tableau 17 : Hiérarchie finale des molécules-cibles sélectionnées (les 26 molécules avant les outputs) Les molécules sont classées en ordre décroissant de la différence arithmétique [sensibilité - spécificité]. On voit que TUP1 est la seule molécule du graphe pour laquelle la sensibilité est supérieure à la spécificité. Ceci l'individualise des autres sur des critères qualitatifs. Ici encore, TUP1 est en première position dans la hiérarchie des cibles potentielles.
La molécule-cible TUP1 est classée en première position. C'est bien la molécule du graphe dont l'inhibition par délétion du gène induit l'évolution du système biologique étudié vers l'état à atteindre. L'efficacité de la méthode est donc vérifiée : sélection d'un nombre restreint de molécules- cibles, et pertinence du classement hiérarchique de ces cibles.
L'ensemble de la mise en œuvre des étapes A et B, ou A', puis C, D et E a été réalisé plus d'une dizaine de fois, et a systématiquement donné des résultats similaires à chaque fois (avec classement de TUP1 toujours dans les 5 premières molécules par le premier mode de classement). Ceci montre sans ambiguïté la reproductibilité de la méthode ainsi que son efficacité.
Si le classement des molécules était effectué au hasard, la probabilité de classer la molécule TUP1 en première position ne serait que de 0,0075 (= 1/133), et la probabilité de répéter 10 fois ce classement de l'ordre de 5.10" 32 (= 0.007510). Il est clair que le résultat obtenu dans cet exemple de mise en œuvre n'est pas dû au hasard et est extrêmement significatif.
Ceci montre la capacité de la méthode à prévoir la cible sur laquelle il faudrait agir, et le type d'action (ici une inhibition) pour atteindre un état donné, c'est à dire la capacité de la méthode à identifier / sélectionner / classer des cibles thérapeutiques.

Claims

REVENDICATIONS
1. Procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, permettant l'analyse dudit réseau d'interactions lorsqu'un stimulus est appliqué au modèle dynamique, en vue notamment de hiérarchiser des molécules biologiques ou de sélectionner des cibles thérapeutiques vis-à-vis d'un problème biologique donné, pour en particulier définir une action thérapeutique à appliquer auxdites molécules, ledit procédé étant mis en œuvre par un système informatique- et comprenant les étapes suivantes : A) à partir d'un graphe statique dont les sommets représentent des molécules biologiques et les arcs représentent des interactions physico-chimiques existant entre ces molécules, associer une variable quantitative Xi mesurée expérimentalement à chaque sommet i, et une relation mathématique à chaque arc du graphe, chacune desdites relations présentant les caractéristiques suivantes : - elle comprend un terme inertiel (i) qui tend vers une limite finie; - elle comprend un terme (ii) tendant à faire revenir les variables Xj à leur état initial, de signe inverse au terme inertiel (i), et dont la variation en fonction du temps croit en valeur absolue de façon plus lente que la variation en fonction du temps du terme inertiel (i); - elle comporte un facteur de pondération w qui permet de tenir compte de la combinaison d'effets pouvant s'exercer sur chaque sommet du graphe; B) calculer les paramètres de chaque relation à partir de données expérimentales quantitatives concernant les sommets du graphe, par la mise en œuvre de techniques d'apprentissage par descente de gradient utilisées pour le paramétrage de réseaux.
2. Procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon la revendication 1 , dans lequel les relations mathématiques associées aux arcs sont continues.
3. Procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon la revendication 1 ou 2, dans lequel chaque variable quantitative Xi associée à un sommet représente la variation relative de la quantité de la molécule correspondant audit sommet, par rapport à la quantité de la même molécule dans un état étalon du système biologique.
4. Procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon l'une des revendications 1 à 3, dans lequel le terme inertiel (i) est exprimé sous la forme d'une relation mathématique présentant une ou plusieurs ênflexion(s).
5. Procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon la revendication 4, dans lequel le terme inertiel (i) est exprimé sous la forme d'une relation sigmoïde ou d'une relation d'oscillation.
6. Procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon l'une quelconque des revendications 1 à 5, dans lequel l'étape B) est effectuée par descente de gradient simple, en prenant comme base de calcul les couples de données (X|, Xj) fournis par les données expérimentales, indépendamment les uns des autres, ou par descente de gradient dans le temps, les couples (Xj, Xj) n'étant alors pas considérés comme indépendants les uns des autres.
7. Procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon l'une quelconque des revendications 3 à 6, dans lequel les données expérimentales quantitatives concernant les sommets du graphe sont obtenues par l'utilisation des techniques de criblage à grande échelle.
8. Procédé d'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique selon l'une quelconque des revendications 3 à 6, dans lequel l'état étalon est un état stable du système biologique, dans lequel la quantité de chaque molécule associée à un sommet du graphe est mesurée expérimentalement.
9. Procédé d'analyse d'un réseau d'interactions moléculaires dans un système biologique, par la mise en œuvre d'un système informatique, comportant les étapes suivantes : A') utilisation d'un modèle dynamique du réseau d'interactions moléculaires, ledit modèle étant construit à partir d'un graphe statique dont les sommets représentent des molécules biologiques du système biologique et les arrêtes représentent des interactions physico- chimiques entre ces molécules, et à partir de données expérimentales concernant les taux ou les activités de ces molécules biologiques et susceptible d'être obtenu, par un procédé selon l'une quelconque des revendications 1 à 8, , C) un état du graphe, mesuré expérimentalement, est choisi comme "état à modifier", et la durée du processus biologique à simuler est définie et découpée en une série de pas de temps, D) plusieurs procédures itératives de simulation sont effectuées, comprenant chacune les étapes suivantes : a) un stimulus est imposé à l'état à modifier, c'est-à-dire que la valeur d'une ou de plusieurs des variables quantitatives associées aux sommets du graphe est modifiée, constituant ainsi un état de départ de la simulation ; b) à partir de l'état de départ de la simulation, un calcul de propagation est effectué au sein du graphe.
10. Procédé de sélection de cibles thérapeutiques mettant en œuvre un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, par la mise en œuvre d'un système informatique, comprenant les étapes et caractéristiques suivantes : A') utilisation d'un modèle dynamique du réseau d'interactions moléculaires, ledit modèle étant construit à partir d!un graphe statique dont les sommets représentent des molécules biologiques du système biologique et les arrêtes représentent des interactions physico- chimiques entre ces molécules et à partir de données expérimentales concernant les taux ou les activités de ces molécules biologiques et susceptible d'être obtenu, par un procédé selon l'une quelconque des revendications 1 à 8, ; C) un état du graphe, mesuré expérimentalement, est choisi comme "état à modifier", et la durée du processus biologique à simuler est définie et découpée en une série de pas de temps, et un état du graphe correspondant à un "état à atteindre" du système biologique est choisi comme "état final du graphe" à atteindre; D) plusieurs procédures itératives de simulation sont effectuées, comprenant chacune les étapes suivantes : a) un stimulus est imposé à l'état à modifier, c'est-à-dire que la valeur d'une ou de plusieurs des variables quantitatives associées aux sommets du graphe est modifiée, constituant ainsi un état de départ de la simulation ; b) à partir de l'état de départ de la simulation, un calcul de propagation est effectué au sein du graphe ; c) un calcul de proximité entre P"état final du graphe " obtenu à l'issue de l'étape b) et l'état à modifier, ou entre P"état final du graphe " et un état voulu est effectué ; E) à partir de l'ensemble des proximités statistiques calculées à l'étape D), les sommets, et les stimuli imposés sur ces sommets, sont hiérarchisés, les sommets hiérarchisés correspondant à des cibles thérapeutiques classées.
11. Procédé selon les revendications 9 ou 10, dans lequel les relations mathématiques associées aux arcs du graphe .à l'étape A') sont continues.
12. Procédé selon l'une quelconque des revendications 9 à 11 , dans lequel le calcul de propagation effectué à l'étape D) b) est effectué pendant un nombre de pas de temps tel que la durée de la simulation n'excède pas la durée du processus biologique à simuler définie à l'étape C).
13. Procédé selon l'une quelconque des revendications 9 à 12, dans lequel les pas de temps définis à l'étape C) sont d'un ordre de grandeur inférieur à celui des durées expérimentales réelles séparant les séries de données expérimentales quantitatives utilisées pour le calcul des paramètres des relations, à l'étape B) d'un procédé selon l'une quelconque des revendications 1 à 8.
14. Procédé de sélection de cibles thérapeutiques selon l'une quelconque des revendications 10 à 13, dans lequel les stimuli imposés à l'étape D) a) concernent, pour chacune des simulations, un sommet unique, et dans lequel le résultat de l'étape E) est un classement des sommets, de celui sur lequel un stimulus est le plus susceptible d'aboutir à l'état voulu à partir de l'état à modifier, jusqu'à celui sur lequel un stimulus est le moins susceptible d'avoir cet effet.
15. Procédé de sélection de cibles thérapeutiques selon l'une quelconque des revendications 10 à 13, comprenant les étapes suivantes : - un premier classement hiérarchique des sommets est obtenu en effectuant les étapes A'), C), D) et E) en imposant, pour chacune des simulations de l'étape D), des stimuli qui concernent un sommet unique ; - une étape supplémentaire D2) est ensuite effectuée, correspondant à l'étape D) dans laquelle les stimuli imposés à chaque simulation sont exercés sur deux sommets, soit en testant toutes les combinaisons de deux sommets possibles, soit en limitant ces calculs aux combinaisons de deux sommets parmi un certain nombre des sommets les mieux classés à l'étape E), - à partir de l'ensemble des proximités statistiques calculées à l'étape D2), une étape supplémentaire E2) de classement hiérarchique des associations de deux sommets sur lesquels des stimuli sont le plus susceptibles d'avoir l'effet voulu, est effectuée.
16. Procédé de sélection de cibles thérapeutiques selon la revendication 15, comportant en outre une étape D3) correspondant à l'étape D) dans laquelle les stimuli imposés à chaque simulation sont exercés sur trois sommets, soit en testant toutes les combinaisons de trois sommets possibles, soit en limitant ces calculs aux combinaisons de trois sommets choisis parmi un certain nombre des sommets les mieux classés à l'étape E) et des combinaisons de deux sommets les mieux classées à l'étape E2), ladite étape D3) étant suivie d'une étape E3) de classement hiérarchique des associations de trois sommets sur lesquels des stimuli sont le plus susceptibles d'avoir l'effet voulu.
17. Procédé de sélection de cibles thérapeutiques selon la revendication 16, comportant en outre une étape D4) correspondant à l'étape D) dans laquelle les stimuli imposés à chaque simulation sont exercés sur quatre sommets, soit en testant toutes les combinaisons de quatre sommets possibles, soit en limitant ces calculs aux combinaisons de quatre sommets choisis parmi un certain nombre des sommets et combinaisons de sommets les mieux classés aux étapes E), E2) et E3), ladite étape D4) étant suivie d'une étape E4) de classement hiérarchique des associations de quatre sommets sur lesquels des stimuli sont le plus susceptibles d'avoir l'effet voulu.
18. Procédé de sélection de cibles thérapeutiques selon la revendication 17, dans lequel les étapes D et E sont répétées de façon itérative en augmentant le nombre de sommets sur lesquels s'exercent les stimuli imposés pour les simulations.
19. Procédé de sélection de cibles thérapeutiques selon l'une quelconque des revendications 10 à 13 et 15 à 18, dans lequel, pour les simulations impliquant des stimuli sur plusieurs sommets, les stimuli sont exercés sur ces différents sommets simultanément ou non.
20. Procédé de sélection de cibles thérapeutiques selon l'une quelconque des revendications 15 à 19, comportant en outre une étape de classement statistique des proximités de graphes de toutes les simulations effectuées, intégrant l'ensemble des classements précédemment obtenus.
21. Procédé selon l'une quelconque des revendications 9 à 20, dans lequel l'étape A') correspond sensiblement aux étapes A) et/ou B) de l'une quelconque des revendications 1 à 8.
22. Procédé selon l'une quelconque des revendications 1 à 21 , dans lequel pour au moins une partie des interactions physico-chimiques entre les molécules du système biologique, la relation entre les variables Xj et Xj, deux à deux est de la forme : w .Xj = mj .(d2Xi / dt2) + 2 Ay .(dXi / dt) + ωy2.Xj, dans laquelle mi .(d2Xi / dt2) + ωy2.Xj correspond au terme inertiel (i), 2 .A .(dXj / dt) correspond au terme de retour à l'état initial (ii), Xi est une variable associée à la molécule i dXj / dt est la dérivée de Xj en fonction du temps d2Xj / dt2 est la dérivée seconde de Xj en fonction du temps Xj est une variable associée à la molécule j, mj représente l'inertie de i, λy régit le retour à l'état d'équilibre de Xi, la pulsation ω correspond au temps de réponse de Xj à la variation de Xj, et W est un facteur de couplage représentant la force de l'interaction entre les molécules i et j, correspondant à une pondération de l'effet de chaque molécule j sur la molécule i vis-à-vis de la résultante de l'ensemble des effets combinés de toutes les molécules j exerçant un effet sur i.
23. Procédé selon l'une quelconque des revendications 1 à 22, dans lequel pour au moins une partie des interactions physico-chimiques entre les molécules du système biologique, la relation entre les variables Xj et Xj, deux à deux est établie par une relation sigmoïde comportant un facteur de retardement associée à une fonction de décroissance linéaire.
24. Procédé selon la revendication 23, dans lequel pour au moins une partie des interactions physico-chimiques entre les molécules du système biologique, la relation entre les variables Xj et Xj, est de la forme : (dXj/dt) = Kn . [ 1 / (1 + e"∑ wij Xj " bi) ] - K2i . Xi , où : le terme sigmoïde Ku . [1 / (1 + e"∑ wlj Xj " b')] correspond au terme inertiel (0, le terme K2j . Xi correspond au terme de retour à l'état initial (ii), avec : X = variable associée au sommet i, Xj = variable associée au sommet j, wij = facteur de couplage représentant la force de l'interaction entre les molécules i et j, correspondant à une pondération de l'effet de chaque molécule j sur la molécule i vis-à-vis de la résultante de l'ensemble des effets combinés de toutes les molécules j exerçant un effet sur i., bi = facteur de retardement, K-ϋ = facteur de limite maximale de variation de Xi, et K2i = facteur de retour à l'équilibre.
25. Procédé selon l'une quelconque des revendications 1 à 24, dans lequel pour au moins une partie des interactions physico-chimiques entre les molécules du système biologique, la relation entre les variables Xi et Xj, est une fonction polynôme de type wy Xj = ∑ bmi.Xim = typ-Di .Xip'1 + ... + b3i .Xi3 + b2i .Xi2 + bu .Xj + boi d'ordre strictement inférieur au nombre p de couples (Xit, Xjt) de valeurs expérimentales du niveau de taux ou d'activité Xj ou Xj des molécules i et j, respectivement, à différents instants t, les paramètres bmj étant calculés à partir des p couples expérimentaux (Xit, Xjt) disponibles, et Wy étant un facteur de couplage représentant la force de l'interaction entre les molécules i et j, correspondant à une pondération de l'effet de chaque molécule j sur la molécule i vis-à-vis de la résultante de l'ensemble des effets combinés de toutes les molécules j exerçant un effet sur i.
26. Procédé selon l'une quelconque des revendications 1 à 25, dans lequel pour au moins une partie des interactions physico-chimiques entre les molécules du système biologique, la relation entre les variables Xi et Xj, est une fonction de type dérivée de polynôme
Figure imgf000127_0001
[m.O^p'-l] avec 1 < p' < p - 1 , p étant le nombre de couples expérimentaux (Xjt, Xjt) disponibles.
27. Procédé selon la revendication 26, dans lequel p'=3.
28. Procédé selon l'une quelconque des revendications 1 à 27, dans lequel pour au moins une partie des molécules du système biologique, la résultante globale de n interactions exercées par des molécules 1 à n sur une molécule i est une somme pondérée des actions des molécules 1 à n sur la molécule i, de la forme
Figure imgf000127_0002
i est la relation associée à l'arc (i, j) pour chaque couple (i, j) et ay = (dXj/dt) / Σ (dXj/dt). D':1 -»n]
29. Procédé selon l'une quelconque des revendications 1 à 27, dans lequel pour au moins une partie des molécules du système biologique, la résultante globale de n interactions exercées par des molécules 1 à n sur une molécule i est une somme pondérée des actions des molécules 1 à n sur la molécule i, de la forme FG (∑ j -» i ) = Σ ay . i, où U:1 -»n] D-1→>n] fji est la relation associée à l'arc (i, j) pour chaque couple (i, j) et ay = (d2Xj/dt2) / Σ (d2Xj/dt2) D-1 -»n]
30. Procédé de détermination du mode d'action d'un xénobiotique, consistant à mettre en œuvre un procédé selon la revendication 9 ou 10 dans les conditions suivantes : (i) le système biologique dans lequel un réseau d'interactions moléculaires est étudié est concerné par l'action du xénobiotique ; (ii) l'"état à modifier" choisi à l'étape C), correspond à un état observé expérimentalement avant l'administration dudit xénobiotique ; (iii) on identifie les modifications à apporter au cours de l'étape D)a) pour que le calcul effectué à l'étape D)b) montre une évolution du système vers un état proche de l'état observé après administration du xénobiotique.
31 . Procédé de prédiction d'éventuels effets indésirables d'un traitement, consistant à mettre en œuvre un procédé selon la revendication 9 ou 10 dans les conditions suivantes : (i) le système biologique dans lequel un réseau d'interactions moléculaires est étudié est concerné par le traitement ; (ii) les modifications de l'étape D)a) correspondent aux modifications des niveaux de taux ou d'activité des molécules cibles observées ou souhaitées lors de l'application du traitement ; (iii) l'étape D)b) de calcul de l'évolution du système biologique est suivie d'une analyse de sous-parties du système correspondant à des fonctions physiologiques connues, afin d'identifier les éventuelles évolutions de ces sous-parties vers des états proches d'états pathologiques de référence.
32. Procédé pour hiérarchiser des cibles thérapeutiques potentielles pour une pathologie, consistant à mettre en œuvre un procédé selon l'une quelconque des revendications 9 à 29 et 31 , puis à déterminer le rapport "bénéfice thérapeutique / effets indésirables" d'une action sur chacune des cibles thérapeutiques potentielles.
33. Procédé selon l'une quelconque des revendications 1 à 32, dans lequel le nombre de variables Xi du réseau d'interactions moléculaires considéré est supérieur à environ 100, supérieur à environ 200, ou supérieur à environ 300.
34. Procédé selon l'une quelconque des revendications 1 à 32, dans lequel le nombre de variables Xj des réseaux d'interactions moléculaires considéré est inférieur à environ 100 et en ce que l'on associe lesdits réseaux d'interactions moléculaires, pour former une association de réseaux.
35. Procédé selon la revendication 34, dans lequel le nombre de réseaux associés est compris entre 2 et environ 100.
36. Utilisation d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique susceptible d'être obtenu par un procédé selon l'une quelconque des revendications précédentes, pour étendre un graphe statique dont les sommets représentent des molécules biologiques et les arcs représentent des interactions physico- chimiques entre ces molécules, pour identifier de nouvelles interactions moléculaires.
37. Système informatique pour l'obtention d'un modèle dynamique d'un réseau d'interactions moléculaires dans un système biologique, et l'analyse de ces interactions moléculaires lorsqu'un stimulus est appliqué au modèle dynamique, comprenant au moins une unité centrale de traitement de données reliée à au moins une base de données expérimentales quantitatives, caractérisé en ce qu'il comprend : A) un module de construction d'un graphe statique, dont les sommets représentent des molécules biologiques et les arcs représentent des interactions physico-chimiques existant entre ces molécules, chaque sommet étant associé à une variable quantitative mesurée expérimentalement et chaque arc du graphe étant associé à une relation mathématique; et C) un module d'apprentissage pour calculer les paramètres de chaque relation à partir des données expérimentales quantitatives concernant les sommets du graphe, par la mise en œuvre de techniques d'apprentissage par descente de gradient utilisées pour le paramétrage de réseaux.
38. Système informatique selon la revendication 37, caractérisé en ce qu'il comprend également : D) un module de simulation pour effectuer plusieurs procédures itératives de simulation consistant à imposer un stimulus à un état de graphe mesuré expérimentalement et choisi comme « état à modifier », le stimulus modifiant la valeur d'une ou de plusieurs des variables quantitatives associées aux sommets du graphe, constituant ainsi un état de départ de la simulation à partir duquel un calcul de propagation est effectué au sein du graphe, pour l'obtention d'un « état final du graphe »; et E) un module d'itération pour la modification du stimulus.
39. Système informatique selon la revendication 38, caractérisé en ce qu'il comprend également : E) un modulé de calcul de proximité entre I' « état final d'un graphe » et I' « état à modifier », ou entre I' « état final d'un graphe » et un état voulu, et de hiérarchisation des sommets et des stimuli imposés sur les sommets du graphe, les sommets hiérarchisés correspondant à des cibles thérapeutiques classées.
PCT/FR2004/002064 2003-08-01 2004-07-30 Methode et systeme de selection de cibles therapeutiques par l'utilisation de reseaux dynamiques d'interactions moleculaires WO2005013173A2 (fr)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CA002534401A CA2534401A1 (fr) 2003-08-01 2004-07-30 Methode et systeme de selection de cibles therapeutiques par l'utilisation de reseaux dynamiques d'interactions moleculaires
EP04786022A EP1649405A2 (fr) 2003-08-01 2004-07-30 METHODE ET SYSTEME DE SELECTION DE CIBLES THERAPEUTIQUES PAR L&apos;UTILISATION DE RESEAUX DYNAMIQUES D INTERACTIONS MOLECULA IRES
US11/342,707 US20060235670A1 (en) 2003-08-01 2006-01-31 Method and system for selecting therapeutic targets using molecular interaction dynamic networks

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR0309557 2003-08-01
FR0309557A FR2858446B1 (fr) 2003-08-01 2003-08-01 Methode d'analyse de reseaux d'interactions moleculaires biologiques

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US11/342,707 Continuation US20060235670A1 (en) 2003-08-01 2006-01-31 Method and system for selecting therapeutic targets using molecular interaction dynamic networks

Publications (2)

Publication Number Publication Date
WO2005013173A2 true WO2005013173A2 (fr) 2005-02-10
WO2005013173A3 WO2005013173A3 (fr) 2005-09-29

Family

ID=34043752

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/FR2004/002064 WO2005013173A2 (fr) 2003-08-01 2004-07-30 Methode et systeme de selection de cibles therapeutiques par l'utilisation de reseaux dynamiques d'interactions moleculaires

Country Status (5)

Country Link
US (1) US20060235670A1 (fr)
EP (1) EP1649405A2 (fr)
CA (1) CA2534401A1 (fr)
FR (1) FR2858446B1 (fr)
WO (1) WO2005013173A2 (fr)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7987234B1 (en) * 2004-10-14 2011-07-26 Oracle America, Inc. Monitoring alert conditions
EP2291791A4 (fr) 2008-05-27 2011-06-22 Sloan Kettering Inst Cancer Modèles de perturbations combinatoires de systèmes biologiques vivants
US10910084B2 (en) 2011-06-01 2021-02-02 Albert-Ludwigs-Universität Freiburg Method for modelling, optimizing, parameterizing, testing and validating a dynamic network with network perturbations
EP2530615A1 (fr) * 2011-06-01 2012-12-05 Albert-Ludwigs-Universität Freiburg Procédé de modélisation, optimisation, paramétrage, test et/ou validation des perturbations de réseau dynamique
US11515004B2 (en) 2015-05-22 2022-11-29 Csts Health Care Inc. Thermodynamic measures on protein-protein interaction networks for cancer therapy
US11456053B1 (en) 2017-07-13 2022-09-27 X Development Llc Biological modeling framework
CN109635439B (zh) * 2018-12-12 2023-05-05 新疆大学 基于传播动力学的扰动环境下生产网络连锁效应研究方法
CN116844645B (zh) * 2023-08-31 2023-11-17 云南师范大学 一种基于多视角分层超图的基因调控网络推断方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002044992A2 (fr) * 2000-11-28 2002-06-06 Physiome Sciences, Inc. Systeme de modelisation de mecanismes biologiques
WO2003017127A2 (fr) * 2001-08-16 2003-02-27 Biotech Research Ventures Pte Limited Modelisation de chemins biochimiques

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002044992A2 (fr) * 2000-11-28 2002-06-06 Physiome Sciences, Inc. Systeme de modelisation de mecanismes biologiques
WO2003017127A2 (fr) * 2001-08-16 2003-02-27 Biotech Research Ventures Pte Limited Modelisation de chemins biochimiques

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
CURTO R ET AL: "Mathematical models of purine metabolism in man" MATHEMATICAL BIOSCIENCES, ELSEVIER, NEW YORK, NY, US, vol. 151, no. 1, juillet 1998 (1998-07), pages 1-49, XP002250783 ISSN: 0025-5564 *
GORYANIN I ET AL: "Mathematical simulation and analysis of cellular metabolism and regulation" BIOINFORMATICS, OXFORD UNIVERSITY PRESS, SURREY, GB, vol. 15, no. 9, septembre 1999 (1999-09), pages 749-758, XP002190906 ISSN: 1367-4803 *
SEDAGHAT AHMAD R ET AL: "A mathematical model of metabolic insulin signaling pathways." AMERICAN JOURNAL OF PHYSIOLOGY, vol. 283, no. 5 Part 1, novembre 2002 (2002-11), pages E1084-E1101, XP002289933 ISSN: 0002-9513 *

Also Published As

Publication number Publication date
FR2858446A1 (fr) 2005-02-04
CA2534401A1 (fr) 2005-02-10
EP1649405A2 (fr) 2006-04-26
US20060235670A1 (en) 2006-10-19
FR2858446B1 (fr) 2007-11-09
WO2005013173A3 (fr) 2005-09-29

Similar Documents

Publication Publication Date Title
Casalino et al. AI-driven multiscale simulations illuminate mechanisms of SARS-CoV-2 spike dynamics
Li et al. Modeling and analysis of RNA‐seq data: a review from a statistical perspective
Hill et al. Inferring causal molecular networks: empirical assessment through a community-based effort
JP6356359B2 (ja) アンサンブルに基づいたリサーチ・レコメンデーションシステムおよび方法
Mathur et al. Gene set analysis methods: a systematic comparison
US10437858B2 (en) Database and data processing system for use with a network-based personal genetics services platform
Ma et al. CGI: a new approach for prioritizing genes by combining gene expression and protein–protein interaction data
Nicolaou et al. De novo drug design using multiobjective evolutionary graphs
Rizopoulos Joint models for longitudinal and time-to-event data: With applications in R
Jaeger et al. Mixed effect models for genetic and areal dependencies in linguistic typology
US8571803B2 (en) Systems and methods for modeling and analyzing networks
Zhang et al. Two heads are better than one: current landscape of integrating QSP and machine learning: an ISoP QSP SIG white paper by the working group on the integration of quantitative systems pharmacology and machine learning
US20110119221A1 (en) Method, system and software arrangement for reconstructing formal descriptive models of processes from functional/modal data using suitable ontology
US20220188657A1 (en) System and method for automated retrosynthesis
CN113012770B (zh) 基于多模态深度神经网络药物-药物相互作用事件预测
Korfmann et al. Deep learning in population genetics
Zhou et al. How do tumor cytogenetics inform cancer treatments? dynamic risk stratification and precision medicine using multi-armed bandits
US20060235670A1 (en) Method and system for selecting therapeutic targets using molecular interaction dynamic networks
McComb et al. Generalized pharmacometric modeling, a novel paradigm for integrating machine learning algorithms: a case study of metabolomic biomarkers
Zhang et al. Minimum conflict individual haplotyping from SNP fragments and related genotype
EP1934873A1 (fr) Procede pour determiner l&#39;etat d&#39;un ensemble de cellules et systeme pour la mise en oeuvre dudit procede
Abe et al. UNMF: a unified nonnegative matrix factorization for multi-dimensional omics data
JP2007505405A (ja) コンピュータモデルを使用して治療標的を同定するための装置および方法
Pe'er From gene expression to molecular pathways
Oloulade et al. Cancer drug response prediction with surrogate modeling-based graph neural architecture search

Legal Events

Date Code Title Description
AK Designated states

Kind code of ref document: A2

Designated state(s): AE AG AL AM AT AU AZ BA BB BG BR BW BY BZ CA CH CN CO CR CU CZ DE DK DM DZ EC EE EG ES FI GB GD GE GH GM HR HU ID IL IN IS JP KE KG KP KR KZ LC LK LR LS LT LU LV MA MD MG MK MN MW MX MZ NA NI NO NZ OM PG PH PL PT RO RU SC SD SE SG SK SL SY TJ TM TN TR TT TZ UA UG US UZ VC VN YU ZA ZM ZW

AL Designated countries for regional patents

Kind code of ref document: A2

Designated state(s): GM KE LS MW MZ NA SD SL SZ TZ UG ZM ZW AM AZ BY KG KZ MD RU TJ TM AT BE BG CH CY CZ DE DK EE ES FI FR GB GR HU IE IT LU MC NL PL PT RO SE SI SK TR BF BJ CF CG CI CM GA GN GQ GW ML MR NE SN TD TG

121 Ep: the epo has been informed by wipo that ep was designated in this application
WWE Wipo information: entry into national phase

Ref document number: 2004786022

Country of ref document: EP

ENP Entry into the national phase

Ref document number: 2534401

Country of ref document: CA

WWE Wipo information: entry into national phase

Ref document number: 11342707

Country of ref document: US

WWP Wipo information: published in national office

Ref document number: 2004786022

Country of ref document: EP

WWP Wipo information: published in national office

Ref document number: 11342707

Country of ref document: US