Kinetics model creation system, metabolism control analysis system, kinetics model creation method, and program

Abstract

The present disclosure proposes a kinetics model generation technique for shortening the total calculation time by increasing the probability that a target parameter value is subjected to parameter value processing while avoiding a situation in which a target parameter value to be collected falls into a local solution. The kinetics model generation technique is a technique for generating a set of target parameter values representing values of unknown parameters that contribute to a kinetics model having a desired property, such generation being on the basis of a non-target parameter value representing the value of an unknown parameter that does not contribute to the kinetics model having the desired property. More specifically, the technique includes: sampling processing for generating a candidate parameter value to be set to an unknown parameter, by generating a random number from a given proposal distribution; rejection determination processing for determining whether or not to accept the candidate parameter value on the basis of the non-target parameter value; and parameter value processing for determining whether or not the kinetics model has the desired property when the candidate parameter value accepted by the rejection determination processing is applied to the kinetics model (see Fig. 3).

Description

Kinetic model creation system, metabolic control analysis system, kinetic model creation method, and program

[0001] The present disclosure relates to a kinetic model creation system, a metabolic control analysis system, a kinetic model creation method, and a program.

[0002] In biological substance synthesis, attempts have been made to improve the productivity of target substances by perturbing specific enzymes in metabolic pathways, such as through mass production via genetic engineering. Metabolic control analysis (MCA) using kinetic models is a well-known method for computationally predicting the effects of perturbations on target substance productivity. Kinetic models describe the dynamics of metabolic pathways specific to each biological species as ordinary differential equations (ODEs), enabling computational simulation of metabolic pathway dynamics through numerical integration. MCA is a theory that quantitatively evaluates the effects of perturbations on the concentrations and fluxes of individual metabolites in metabolic pathways and can be performed computationally on kinetic models. However, predicting the effects of perturbations using MCA requires that the kinetic model possess various properties. For example, to perform MCA, the kinetic model must be able to achieve an asymptotically stable steady state (hereinafter referred to as stable). Furthermore, a high degree of reproducibility of actual cell data by the kinetic model is required as an indicator of the reliability of the prediction of the effects of perturbations. The nature of these kinetic models is V ｍａｘ Ya K Ｍ Since the kinetics depends on the values ​​of parameters such as α, β ...

[0003] Non-Patent Document 1 discloses that a technique called ensemble modeling (EM) is used to collect a large number of target parameter values ​​in a kinetic model having unknown parameters, and to create a set of kinetic models in which the collected parameter values ​​are set as unknown parameters. Note that since collecting parameter values ​​and creating a kinetic model in which the collected parameter values ​​are set as unknown parameters are essentially synonymous, in the following description, collecting parameter values ​​and creating a kinetic model may be interpreted interchangeably.

[0004] 1 is a flowchart showing an outline of parameter value collection using EM. First, in a sampling step 11, random numbers are generated as parameter values ​​from a probability distribution called a proposal distribution (parameter values ​​are sampled). Next, in a parameter value processing step 12, the following processes (i) to (iii) are performed.

[0005] (i) The parameter values ​​sampled in the sampling step 11 are set as unknown parameters of the kinetic model. (ii) It is determined whether the kinetic model of (i) has the desired properties. In other words, it is determined whether the set parameter values ​​are the desired parameter values. (iii) If the set parameter values ​​are the desired parameter values ​​as a result of (ii), they are recorded.

[0006] Next, in the termination determination step 13, if the number of recorded target parameter values ​​reaches a predetermined target number, the series of processes is terminated, and if the target number is not reached, the process returns to the sampling step 11.

[0007] The above series of steps from the sampling step 11 to the termination determination step 13 are repeated many times to collect the target number of target parameter values. The number of repetitions can reach 1,000,000 times. As an example of analysis using a set of collected parameter values, there is a method in which MCA is performed on a set of kinetic models for which each parameter value is set, and the results obtained for each model are statistically analyzed, and this method is also implemented in Non-Patent Document 1.

[0008] Non-Patent Document 2 discloses a method called REKINDLE for generating stable parameter values ​​using a generative deep learning model. This method trains a generative deep learning model that generates parameter values ​​using training data of target parameter values ​​and non-target parameter values ​​collected in the same manner as EM. Figure 2 is a flowchart showing an outline of kinetic model creation using REKINDLE. The sampling step 11, parameter value processing step 12, and termination determination step 13 are the same as those in EM. However, in the parameter value processing step 12, the target parameter values ​​and non-target parameter values ​​are recorded as training data 21. Furthermore, once a target number of training data has been collected, the termination determination step 13 transitions to a training step 22. In the training step 22, the generative deep learning model is trained using the training data 21. In the generative task execution step 23, a large number of target parameter values ​​are generated in a short period of time using a generative deep learning model that has become capable of generating parameter values ​​with a high probability of being the target parameter value through training in the training step 22.

[0009] Nishiguchi, H., Hiasa, N., Uebayashi, K., Liao, J., Shimizu, H., & Matsuda, F. (2019). Transomics data-driven, ensemble kinetic modeling for system-level understanding and engineering of the cyanobacteria central metabolism. Metabolic engineering, 52, 273-283. Choudhury, S., Moret, M., Salvy, P., Weilandt, D., Hatzimanikatis, V., & Miskovic, L. (2022). Reconstructing kinetic models for dynamical studies of metabolism using generative adversarial networks. Nature Machine Intelligence, 4(8), 710-719.

[0010] In EM, the parameter value processing step 12 often involves numerical integration, which requires a long calculation time. On the other hand, the probability of obtaining the target parameter value as a result of the determination in the parameter value processing step 12 is generally low. Therefore, in order to obtain the target number of target parameter values, the parameter value processing step 12 must be performed a number of times that far exceeds the target number. For these reasons, EM has a problem in that when the parameter value processing step 12 involves numerical integration, the total calculation time is very long, sometimes reaching several weeks.

[0011] Regarding this issue, if a small number of target parameter values ​​are identified in the EM process and then the surrounding values ​​are locally sampled, the probability that the target parameter values ​​will be provided to the parameter value processing step 12 may be improved, potentially solving the problem of long EM calculation times. However, this means that the collected parameter values ​​will fall into a local solution, increasing the possibility that parameter values ​​that more closely reproduce the actual metabolic pathway will not be collected.

[0012] Furthermore, as shown in Non-Patent Document 2, REKINDLE makes it possible to collect a large number of target parameter values ​​in a short time in the generation task execution step 23. However, because the collection of training data 21 is performed in the same manner as EM, the same calculation time as EM is already required at the stage before the training step 22. Therefore, it has the problem of requiring a longer calculation time than EM.

[0013] In view of such circumstances, the present disclosure proposes a technique for reducing the total calculation time by increasing the probability that the collected target parameter values ​​are used in parameter value processing while avoiding the target parameter values ​​from falling into a local solution.

[0014] In order to solve the above problems, the present disclosure proposes a kinetic model generation system that collects target parameter values ​​for a kinetic model having unknown parameters and generates a kinetic model by setting the collected target parameter values ​​to the values ​​of the unknown parameters, the kinetic model comprising: a storage device that stores a program for generating a set of target parameter values ​​representing values ​​of unknown parameters that contribute to the kinetic model having desired properties, based on untarget parameter values ​​representing values ​​of unknown parameters that do not contribute to the kinetic model having desired properties; and a processor that reads and executes the program from the storage device, wherein the processor performs: a sampling process that generates candidate parameter values ​​to be set for the unknown parameters by generating random numbers from a given proposal distribution; a rejection determination process that determines whether or not to accept the candidate parameter values, based on the untarget parameter values; and a parameter value process that determines whether or not the kinetic model has desired properties when the candidate parameter values ​​accepted by the rejection determination process are applied to the kinetic model.

[0015] Further features related to the present disclosure will become apparent from the description of this specification and the accompanying drawings. Also, aspects of the present disclosure are achieved and realized by the elements and combinations of various elements and the aspects of the following detailed description and the appended claims. The description of this specification is merely exemplary and does not limit the scope or application of the claims of the present disclosure in any way.

[0016] According to the technology disclosed herein, it is possible to prevent the collected target parameter values ​​from falling into a local solution, increase the probability that the target parameter values ​​are used in parameter value processing, and shorten the total calculation time required to find parameter values ​​that are suitable for the purpose.

[0017] 10 is a flowchart showing an overview of kinetic model creation using EM. FIG. 11 is a flowchart showing an overview of kinetic model creation using REKINDLE. FIG. 12 is a flowchart showing an overview of kinetic model creation processing according to this embodiment. FIG. 13 is a diagram showing an example of the hardware configuration of a kinetic model creation and metabolic simulation system 40 according to this embodiment. FIG. 14 is a diagram showing an example of the software configuration of a kinetic model creation and metabolic simulation system 40 according to this embodiment. FIG. 15 is a flowchart for explaining details of kinetic model creation and metabolic simulation according to this embodiment. FIG. 16 is a diagram showing an example of the configuration of a GUI when input information 001 is received by the input receiving unit 100 according to this embodiment. FIG. 17 is a diagram showing an example of the configuration of a GUI when a list 403 of stable parameter values ​​is output by the parameter value output unit 402 according to this embodiment. FIG. 18 is a diagram showing an example of the configuration of a GUI when a perturbation effect prediction result 505 is output by the prediction result output unit 504 according to this embodiment. FIG. 19 is a flowchart for explaining an example algorithm 1 of a rejection judgment function according to this embodiment. FIG. 19 is a diagram schematically showing the concept of rejection judgment by the processes of S10002 to S10006 of FIG. 10. FIG. 19 is a flowchart for explaining an example algorithm 2 of a rejection judgment function according to this embodiment. 17 is a diagram schematically showing the concept of rejection determination by the processing of S12002 to S12005 in FIG. 12. FIG. 18 is a flowchart for explaining algorithm example 3 of the rejection determination function according to this embodiment. FIG. 19 is a diagram schematically showing the concept of rejection determination by the processing of S14002 to S14004 in FIG. 14. FIG. 19 is a graph showing how stable parameter values ​​are collected over the course of calculation time according to an example. FIG. 20 is a graph showing how unstable parameter values ​​are recorded in the process of collecting stable parameter values ​​in FIG. 16 according to an example. FIG. 21 is a graph showing the final distribution of parameter values ​​determined to be stable and unstable in FIG. 16 and FIG. 17 according to an example.

[0018] In this embodiment, a new technique called empirical rejection is introduced between the sampling step 11 and the parameter value processing step 12 in EM, in which the parameter value processing step 12 involves numerical integration. This technique prevents the collected target parameter values ​​from falling into a local solution, while increasing the probability that the parameter values ​​provided to the parameter value processing step 12 are actually the target parameter values. This proposes a kinetic model creation technology that can shorten the total calculation time. Hereinafter, this embodiment will be described with reference to the accompanying drawings. Note that in the accompanying drawings, functionally identical elements may be denoted by the same numbers. Note that the accompanying drawings show specific embodiments and implementation examples in accordance with the principles of the present disclosure, but these are intended to aid in understanding the present disclosure and should in no way be used to limit the interpretation of the present disclosure.

[0019] Furthermore, although the present embodiment has been described in sufficient detail to enable those skilled in the art to implement the present disclosure, it should be understood that other implementations and forms are possible, and that changes to the configuration and structure and substitutions of various elements are possible without departing from the scope and spirit of the technical ideas of the present disclosure. Therefore, the following description should not be interpreted as being limited thereto.

[0020] <Definitions of terms used in this specification> In this specification, "enzyme" refers to a protein that catalyzes a chemical reaction. "Metabolic pathway" refers to a chain of enzymatic reactions that convert compounds one after another. "Metabolite" is a general term for the starting material, intermediate product, and final product of a metabolic pathway. "Flux" refers to the net rate obtained by subtracting the rate of the forward reaction from the rate of the reverse reaction in a chemical reaction between two types of metabolite. "Productivity" refers to the magnitude of the flux that produces a certain substance.

[0021] "Perturbation" generally means disturbance or change, and here refers to artificially changing the amount or activity of metabolic enzymes. "Rate equation" refers to an equation that describes the rate of a reaction in a metabolic pathway as a function of the concentration of each substance in the metabolic pathway. "Parameter" here refers to the V in the rate equation. ｍａｘ (maximum reaction rate) and K ＭA "parameter value" refers to a constant that determines the properties of a metabolic pathway, such as the Michaelis constant. Here, "parameter value" refers to the value of a single parameter, and to a combination of these values ​​when multiple parameters are considered. "Metabolic Control Analysis (MCA)" refers to a theory for quantitatively evaluating the effect of parameter values, such as enzyme amounts, on each flux and each metabolic substance concentration in a metabolic pathway. In MCA, "flux control coefficient (FCC)" refers to an index that represents the degree to which a certain flux changes with a change in a certain parameter value. Here, "kinetic model" refers to an ODE constructed by combining rate equations, with the concentration of each substance in a metabolic pathway as the dependent variable and time as the independent variable. Here, "having a desired property" of a kinetic model means that the kinetic model has the desired property.

[0022] A "target parameter value" refers to a parameter value set as an unknown parameter that causes the kinetic model to have the desired properties (i.e., contributes to the kinetic model having the desired properties). An "untargeted parameter value" refers to a parameter value set as an unknown parameter that is not a target parameter value (i.e., does not contribute to the kinetic model having the desired properties). A "state" here refers to a combination of values ​​of dependent variables in a kinetic model. The "relative rate of change" of a dependent variable in a kinetic model here refers to the value obtained by dividing the first-order differential coefficient at a certain state by the value at that state. A "steady state" here refers to a state in a kinetic model where the relative rates of change of all dependent variables are zero. The "distance" between one state and another state refers to the distance when the states are viewed as points in space. A steady state in a kinetic model is said to be "asymptotically stable" if, as time progresses, any state that is a finite distance away from that steady state always reaches that steady state. "Linear stability analysis" of a steady state of a kinetic model refers to an analysis that determines that the steady state of the kinetic model is asymptotically stable when the real parts of all eigenvalues ​​of the Jacobian matrix are negative.

[0023] Here, a kinetic model is "stable" if it can be confirmed that the kinetic model has a steady state and that the steady state is asymptotically stable. Here, a kinetic model is "unstable" if it is not stable. A "stable parameter value" refers to a parameter value that makes the kinetic model stable among the parameter values ​​set for the unknown parameters. An "unstable parameter value" refers to a parameter value that is not a stable parameter value.

[0024] Here, "proposal distribution" refers to the probability distribution that a random number follows when it is generated as a parameter value. "Sampling" refers to generating a random number as a parameter value from the proposal distribution. "Sampling function" refers to a function that, when executed, returns a random number that follows the proposal distribution.

[0025] "Rejection" here means that in a series of processing flows including sampling, the sampled random number is not used for subsequent processing, and sampling is performed again. "Acceptance" here means that in a series of processing flows including sampling, the random number whose parameter value has been sampled is not rejected, and is used for subsequent processing. "Rejection determination function" here means a function that takes the sampled parameter value and an unintended parameter value as arguments and returns "reject" or "accept."

[0026] <Outline of Kinetic Model Creation> Figure 3 is a flowchart showing an outline of kinetic model creation according to this embodiment. While kinetic model creation according to this embodiment includes all of the steps included in conventional EM, it is characterized in that a rejection determination step 31 is introduced between the sampling step 11 and the parameter value processing step 12, and that unintended parameter values ​​32 identified in the parameter value processing step 12 are recorded in a storage device 42 (see Figure 4). The rejection determination step 31 determines whether to redo the sampling step 11 or proceed to the parameter value processing step 12 based on the unintended parameter values ​​32 recorded in the parameter value processing step 12. This new technique for determining whether to reject or accept parameter values ​​based on already identified unintended parameter values ​​is called the empirical rejection method. In other words, kinetic model creation according to this embodiment introduces the empirical rejection method into conventional EM.

[0027] In this embodiment, a system is described in which stable parameter values ​​are collected as target parameter values ​​for kinetic models having unknown parameters, and MCA is performed using a set of kinetic models set with the collected stable parameter values ​​to predict the effect of perturbation.

[0028] <Example of Hardware Configuration of Kinetic Model Creation and Metabolic Simulation System> FIG. 4 is a diagram showing an outline of a kinetic model creation and metabolic simulation system 40 according to this embodiment.

[0029] The kinetic model creation and metabolic simulation system 40 can be configured by a computer system, and includes, for example, a control device 41, a storage device 42, an input device 43, and an output device 44.

[0030] The control device 41 reads the kinetic model creation and metabolic simulation program (see Figure 5) from the storage device 42 and expands it into internal memory (not shown), thereby realizing an input receiving unit 100, a parameter value generating unit 200, a parameter value processing unit 300, a termination control unit 400, and a perturbation effect predicting unit 500.

[0031] The storage device 42 stores a kinetic modeling and metabolic simulation program (see FIG. 5) that executes a process of predicting the effect of perturbations on metabolic pathways using an empirical rejection method (see FIG. 3) and a set of collected parameter values. The kinetic modeling and metabolic simulation program corresponds to the flowchart in FIG. 6, which will be described later.

[0032] The input device 43 is composed of, for example, a keyboard, a mouse, a switch, etc., and can be used, for example, by an operator (user) to input information such as a kinetic model having unknown parameters, a target set size N, a proposal distribution, and a rejection judgment function, as well as various instructions, into the kinetic model creation and metabolic simulation system 40.

[0033] The output device 44 is composed of a display device, a printer, etc., and displays, for example, the results of kinetic model creation and metabolic simulation (such as a set of collected parameter values ​​and predicted effects of perturbations, as described below) on a display screen or prints them on a substrate.

[0034] In FIG. 4, the computer system is configured such that the control device 41, the storage device 42, the input device 43, and the output device 44 are connected via data lines / control lines, but the kinetic model creation and metabolic simulation system 40 may also be configured such that each device is remotely located and connected via a communication network.

[0035] <Example of Software Configuration of Kinetic Model Creation and Metabolic Simulation System> Fig. 5 is a diagram showing an outline of the software configuration of a kinetic model creation and metabolic simulation system 40 according to this embodiment. The kinetic model creation and metabolic simulation system 40 mainly comprises, as its software configuration, an input receiving unit 100, a parameter value generating unit 200, a parameter value processing unit 300, a termination control unit 400, and a perturbation effect predicting unit 500. Each of these will be described below.

[0036] (i) Input Receiving Unit 100 The input receiving unit 100 receives, as input information 001, a definition of a kinetic model having unknown parameters, a target set size N, a definition of a proposal distribution for generating random numbers as parameter values, and a definition of a rejection judgment function.

[0037] Next, the input receiving unit 100 generates a sampling function from the definition of the received proposal distribution for sampling in the sampling unit 201. Furthermore, the input receiving unit 100 generates a rejection determination function from the definition of the received rejection determination function for making a rejection determination in the rejection determination unit 202.

[0038] (ii) Parameter Value Generation Unit 200 The parameter value generation unit 200 includes, as its internal configuration, a sampling unit 201 and a rejection determination unit 202. The parameter value generation unit 200 performs processing to generate parameter values ​​to be processed in the parameter value processing unit 300. The sampling unit 201 executes the sampling function generated by the input reception unit 100, and generates random numbers as parameter values.

[0039] Furthermore, rejection determination unit 202 inputs the parameter values ​​generated by sampling unit 201 and list 306 of unstable parameter values ​​to the rejection determination function generated by input receiving unit 100, and obtains a determination result of "reject" or "accept." If the determination result is "reject," parameter value generation processing in sampling unit 201 is performed again, and if the determination result is "accept," processing by parameter setting unit 301 is executed.

[0040] (iii) Parameter Value Processing Unit 300 The parameter value processing unit 300 includes, as its internal configuration, a parameter setting unit 301, a numerical integration unit 302, a steady state determination unit 303, a linear stability analysis unit 304, and a stability recording unit 305.

[0041] The parameter setting unit 301 sets the parameter values ​​(accepted values) generated by the parameter value generation unit 200 to the unknown parameters in the definition of the kinetic model having unknown parameters accepted by the input accepting unit 100, and creates a kinetic model (ordinary differential equation).

[0042] The numerical integration unit 302 performs numerical integration on the kinetic model created by the parameter setting unit 301, and acquires the state after the final step (see the definitions of terms, etc., above). Note that the initial value and the integration interval of the numerical integration may be included in the input information 001 by the operator.

[0043] The steady state determination unit 303 calculates the absolute value of the relative rate of change of each metabolite concentration in that state from the kinetic model and the state after the final step obtained by the numerical integration unit 302, and determines that the state is steady if all of these values ​​are below a predetermined threshold and can be considered close to zero, otherwise determines that the state is not steady. If it is determined that the state is not steady, the processing by the linear stability analysis unit 304 is skipped and processing by the stability recording unit 305 is executed. The threshold value for determination may be included in the input information 001.

[0044] The linear stability analysis unit 304 performs a linear stability analysis on the kinetic model in the steady state obtained by the numerical integration unit 302, and determines whether the steady state is asymptotically stable. Alternatively, instead of verifying that the real part of the eigenvalue is negative in the linear stability analysis, the linear stability analysis unit 304 may determine whether the steady state can be considered asymptotically stable by verifying that the real part is below a predetermined threshold close to zero.

[0045] The stability recording unit 305 adds the parameter value generated by the parameter value generating unit 200 to the list of stable parameter values ​​403 when the determination result of the steady state determining unit 303 is "considered to be in a steady state" and the determination result of the linear stability analyzing unit 304 is "asymptotically stable." On the other hand, if this is not the case, the stability recording unit 305 adds the parameter value to the list of unstable parameter values ​​306. The list of unstable parameter values ​​306 can be stored in an internal memory (not shown) or the storage device 42. Note that unstable parameter values ​​are basically acquired and accumulated in the parameter value processing step 12, but parameter values ​​that are clearly not intended may be included in the list of unstable parameter values ​​306 in advance.

[0046] In this embodiment, the target property (desired property) of the kinetic model is "stability," but the target property is not limited to "stability." As an example, the target property may be determined to be a stable kinetic model and the obtained steady state being able to reproduce actual cell data to a certain degree or higher. In this case, an example of a method for determining whether the kinetic model has the target property is to, if the steady state determination unit 303 and the linear stability analysis unit 304 determine that the kinetic model is stable, subsequently calculate the value of an evaluation function (a function that evaluates the degree to which the simulation results and actual data match), such as mean square error, root mean square error, cross entropy error, mean absolute error, or squared logarithmic error, using the obtained steady state and actual cell data. If the value is below a predetermined threshold, the target property is determined to be present; otherwise, the target property is determined to be absent. As another example, the target property may be determined to be a state transition obtained as a result of numerical integration being able to reproduce actual cell time-series data to a certain degree or higher. In this case, an example of a method for determining whether a kinetic model has the desired properties is a method in which, instead of processing by the steady state determination unit 303 and the linear stability analysis unit 304 in this embodiment, the value of the aforementioned evaluation function is calculated using the time series of the state obtained by numerical integration in the numerical integration unit 302 and the time series data of the actual cell, and if the value is below a predetermined threshold, it is determined that the kinetic model has the desired properties, and if not, it is determined that the kinetic model does not have the desired properties.

[0047] (iv) Termination Control Unit 400 The termination control unit 400 includes, as its internal components, a termination determination unit 401 and a parameter value output unit 402. The termination determination unit 401 determines whether the number of parameter values ​​recorded in the stable parameter value list 403 has reached the target set size N received by the input receiving unit 100. If the number of parameter values ​​has not reached the target set size N, the process returns to the processing by the sampling unit 201. If the number of parameter values ​​has reached the target set size N, the parameter value output unit 402 executes processing.

[0048] The parameter value output unit 402 outputs a list of stable parameter values ​​403. The list of stable parameter values ​​403 can be stored in a stable parameter value DB (not shown), displayed on a display screen, or output by a printer, for example.

[0049] (v) Perturbation Effect Prediction Unit 500 The perturbation effect prediction unit 500 executes perturbation effect prediction processing, and includes, as its internal configuration, a parameter setting unit 501, an MCA unit 502, a statistical processing unit 503, and a prediction result output unit 504.

[0050] The parameter setting unit 501 sets each of the parameter values ​​included in the stable parameter value list 403 as an unknown parameter in the definition of the kinetic model having unknown parameters received by the input receiving unit 100 .

[0051] The MCA unit 502 performs MCA on the kinetic model for which unknown parameters have been set by the parameter setting unit 501, and calculates the FCC for each combination of the perturbation target and the reaction affected by the perturbation.

[0052] The statistical processing unit 503 calculates statistics by aggregating the FCCs calculated by the MCA unit 502 among the parameter values ​​used. The statistics include, for example, the median, the mean, the variance, or a combination thereof.

[0053] The prediction result output unit 504 outputs the FCC statistics calculated by the statistical processing unit 503 for each combination of perturbation target and flux as perturbation effect prediction results 505. The perturbation effect prediction results 505 may be sorted based on the magnitude of the FCC statistics, and may be stored in a perturbation effect prediction result DB (not shown), displayed on a display screen, or output by a printer, for example.

[0054] <Details of Kinetic Model Creation and Metabolic Simulation> Fig. 6 is a flowchart for explaining the details of kinetic model creation and metabolic simulation according to this embodiment. In the explanation of each step described below, the operating subject is each processing unit (input receiving unit 100, sampling unit 201, etc.) shown in Fig. 5, but since each processing unit is realized by a program deployed in control device 41, the operating subject of each step may also be read as control device 41.

[0055] (i) Step S6001: The input receiving unit 100 receives, as input information 001, a definition of a kinetic model having unknown parameters, a target set size N, a definition of a proposal distribution for generating random numbers as parameter values, and a definition of a rejection judgment function. The input information 001 is, for example, information input by an operator via the input device 43. Here, the definition of the kinetic model having unknown parameters is input, for example, as a file written in SBML (Systems Biology Markup Language) format. Furthermore, the definition of the proposal distribution is input, for example, as a source file describing a function for generating random numbers based on the proposal distribution, which is then read so that it can be executed in subsequent processing, thereby generating a sampling function. The definition of the rejection determination function is input, for example, as a source file (a file describing the contents of the rejection determination process: the user may specify the rejection criteria, or the file may be built into the system) describing a function that takes the parameter value to be determined and the list 306 of unstable parameter values ​​as arguments and returns a determination result of "reject" or "accept", and is read so that it can be executed in later processing, generating the rejection determination function (described in detail below: see FIG. 10 etc.). A natural number is input as the target set size N.

[0056] (ii) Step S6002 The sampling unit 201 executes the sampling function generated in step S6001 to generate random numbers as parameter values ​​(candidate parameter values).

[0057] (iii) Step S6003 The rejection determination unit 202 inputs the parameter value generated in step S6002 and the parameter value included in the list 306 of unstable parameter values ​​into the rejection determination function, and obtains a determination result of "reject" or "accept."

[0058] (iv) Step S6004: The rejection determination unit 202 determines whether the determination result in step S6003 is "rejection." If the determination result is "rejection" (Yes in step S6004), the process returns to step S6002. If the determination result is "acceptance" (No in step S6004), the process proceeds to step S6005.

[0059] (v) Step S6005 The parameter setting unit 301 sets the parameter value determined as "accepted" in step S6004 to the unknown parameter in the definition of the kinetic model having the unknown parameter accepted in step S6001, and creates a kinetic model.

[0060] (vi) Step S6006: The numerical integration unit 302 performs numerical integration on the kinetic model created in step S6005 to obtain the state after the final step. The initial value and integration interval of the numerical integration may be included in the input information 001. Note that, since the kinetic model is a mathematical model expressed by a differential equation, when the state at each time is input, the degree of change (progression) in the concentration of each metabolite (rate of change) is output.

[0061] The steady state determination unit 303 calculates the absolute value of the relative rate of change of each metabolite concentration in that state (because the relative rate of change may be negative) from the kinetic model and the state after the final step obtained by the numerical integration unit 302. The steady state determination unit 303 then determines that the state is steady if all of the calculated relative rates of change are below a predetermined threshold and can be considered close to zero; otherwise, it does not determine that the state is steady. If it is determined that the state is not steady, the processing by the linear stability analysis unit 304 is skipped, and the process proceeds to step S6007. The threshold for determination may be included in the input information 001.

[0062] The linear stability analysis unit 304 performs a linear stability analysis of the kinetic model in the steady state obtained by the numerical integration unit 302, and determines whether the steady state is asymptotically stable. Alternatively, instead of verifying that the real part of the eigenvalue is negative in the linear stability analysis, the linear stability analysis unit 304 may determine whether the steady state can be considered asymptotically stable by verifying that the real part is below a predetermined threshold close to zero.

[0063] (vii) Step S6007 The stability recording unit 305 adds the parameter value determined as "accepted" in step S6004 to the list of stable parameter values ​​403 if the determination result of the steady state determination unit 303 in step S6006 is "considered to be in a steady state" and the determination result of the linear stability analysis unit 304 is "asymptotically stable." On the other hand, if the determination result is not true, the stability recording unit 305 adds the parameter value to the list of unstable parameter values ​​306. The list of unstable parameter values ​​306 is stored in an internal memory (not shown), the storage device 42, etc.

[0064] (viii) Step S6008 The termination determination unit 401 determines whether the number of parameter values ​​recorded in the list of stable parameter values ​​has reached the target set size N accepted by the input acceptance unit 100. If the number of parameter values ​​has reached the target set size N (Yes in step S6008), the process proceeds to step S6009. If the number of parameter values ​​has not reached the target set size N (No in step S6008), the process returns to step S6002.

[0065] The parameter value output unit 402 outputs a stable parameter value list 403 in which N parameter values ​​are recorded. The stable parameter value list 403 can be stored in a stable parameter value DB (not shown), displayed on a display screen, or output by a printer, for example.

[0066] (ix) Step S6009 The parameter setting unit 501 sets one of the parameter values ​​included in the list 403 of stable parameter values ​​as the unknown parameter in the definition of the kinetic model having the unknown parameter received in step S6001.

[0067] (x) Step S6010 The MCA unit 502 performs MCA on the kinetic model in which the unknown parameters have been set in step S6009, and calculates the FCC for each combination of the perturbation target and the reaction affected by the perturbation.

[0068] As described above, by the processing from step S6009 to step S6010, an FCC is calculated for each combination of the perturbation target and the reaction affected by the perturbation for each parameter value. By performing this for each of the N parameter values ​​included in the list of stable parameter values, N FCCs are calculated for each combination of the perturbation target and the reaction affected by the perturbation. Once the FCCs have been calculated for all N parameter values, the processing proceeds to step S6011. (xi) Step S6011

[0069] The statistical processing unit 503 calculates a statistic by aggregating the FCCs calculated for the N parameter values ​​among the parameter values ​​(the statistic can be regarded as a representative value of the FCCs for the N parameters). In this embodiment, the median is calculated as the statistic, but it may also be, for example, the mean, the variance, or a combination thereof.

[0070] The prediction result output unit 504 outputs the FCC statistics calculated by the statistical processing unit 503 for each combination of perturbation target and flux as perturbation effect prediction results 505. The perturbation effect prediction results 505 may be sorted based on the magnitude of the FCC statistics, and may be stored in a perturbation effect prediction result DB (not shown), displayed on a display screen, or output by a printer, for example.

[0071] 7 to 9 are diagrams showing examples of GUI (Graphical User Interface) display operated by an operator in the kinetic model creation and metabolic simulation system 40 according to this embodiment. Fig. 7 is a diagram showing an example of GUI display when input information 001 is received by the input receiving unit 100. Fig. 8 is a diagram showing an example of display when a list 403 of stable parameter values ​​is output by the parameter value output unit 402. Fig. 9 is a diagram showing an example of display when a perturbation effect prediction result 505 is output by the prediction result output unit 504.

[0072] (i) As shown in FIG. 7, when the operator refers to a file, inputs a natural number as input information 001 on the input screen 700, and presses the [Execute] button, the control device 41 starts the processing consisting of the parameter value generation unit 200, the parameter value processing unit 300, and the termination control unit 400, as described above.

[0073] (ii) When execution is instructed via the input screen 700 and the initiated processing is completed, the control device 41 displays a list 403 of stable parameter values ​​on a stable parameter value list output screen 800, as shown in FIG. 8 . Note that FIG. 8 illustrates a case where there are multiple dimensions (types of unknown parameters) of parameter values, with each row in the table representing one (one set) of parameter values. Furthermore, when a specific reaction is selected in the illustrated metabolic pathway (including multiple reactions), the rate equation for the selected reaction is displayed, and the parameters included in the displayed rate equation are highlighted in the table (e.g., Vmax1 and KM1 in FIG. 8 ). Furthermore, by pressing the [Save] button, the operator can save the displayed list 403 of stable parameter values ​​as a file in a format such as CSV (Comma-Separated Values). Furthermore, when the operator presses the [Perturbation Effect Prediction] button, the control device 41 starts processing performed by the perturbation effect prediction unit 500.

[0074] Upon completion of the process initiated by pressing the [Perturbation Effect Prediction] button on the stable parameter value list output screen 800, the control device 41 displays the perturbation effect prediction results 505 in table format on the perturbation effect prediction result output screen 900, as shown in FIG. 9 . In FIG. 9 , the perturbation effect prediction unit 500 calculates the median FCC between parameter values, and the perturbation target enzymes are sorted and displayed according to the magnitude of the median. For example, by referencing the column of reactions affected by perturbation in the table, which column contains the reaction (target reaction) for which the flux is desired to be increased, the operator can identify the perturbation target with the largest median FCC, i.e., the enzyme predicted to increase the flux of the target reaction when increased. Furthermore, the operator can save the displayed perturbation effect prediction results 505 as a file in a format such as CSV by pressing the [Save] button.

[0075] <Example of Algorithm of Rejection Function> An example of the algorithm of the rejection determination function included in the input information 001 will now be described.

[0076] (1) Algorithm Example 1: Rejection Decision Algorithm Using Probability Density Function FIG. 10 is a flowchart for explaining algorithm example 1 of the rejection decision function.

[0077] (i) Step S10001 The rejection determination unit 202 of the parameter value generation unit 200 receives as input a parameter value q to be determined and a list Q of parameter values ​​that are not of interest. ｕｎ Then, the rejection determination unit 202 receives the list of unstable parameter values ​​Q ｕｎ Take one parameter value from q ｕｎ，ｉ Then, the following processing from step S10002 to step S10005 is performed.

[0078] (ii) Step S10002 The rejection determination unit 202 calculates the average q ｕｎ，ｉ , variance σ ２ The probability density function of the normal distribution f ｉ (x), where σ ２ The specific value of is part of the rejection function and is not received as an input, but σ ２The kinetic model creation and metabolic simulation system 40 and the rejection judgment function may be configured and defined so that the value of (i) can be included in the input information 001 and input to the rejection judgment function. In the case of a probability distribution with multiple parameters, at least one of the parameters (e.g., mean or variance) may be set to an unintended parameter value, while the remaining parameters may be set to values ​​input from an external source. For example, of the mean and variance parameters of the normal distribution, the mean may be set to an unintended parameter value, while the variance may be set to a value input from an external source. Furthermore, the value does not have to be input from an external source; in that case, it can be included as a constant in the definition of the rejection judgment function. Note that while a normal distribution is used here as an example of a probability distribution, a log-normal distribution, a uniform distribution, a gamma distribution, a triangular distribution, or a mixture of these may also be used instead of the normal distribution.

[0079] (iii) Step S10003 The rejection determination unit 202 ｉ = f ｉ (q) / f ｉ (q ｕｎ，ｉ ) is calculated.

[0080] (iv) S10004 The rejection determination unit 202 generates a random number X that follows a continuous uniform distribution U(0,1). ｉ Generates.

[0081] (v) S10005 The rejection determination unit 202 ｉ >X ｉ It is determined whether or not p ｉ >X ｉ If pi≦Xi (if No in step S10005), the rejection determination unit 202 determines Q ｕｎ to another parameter value q ｕｎ，ｉ＋１ is extracted and the process returns to step S10002.

[0082] As described above, by the processing from step S10002 to step S10005, q ｕｎ，ｉ This series of processes is performed using Q ｕｎIf the process proceeds to step S10006 due to the conditional branch in step S10005 during the repetition, the parameter values ​​that have not yet been used in the process are ｕｎ If the process never proceeds to step S10006 due to the conditional branch in step S10005 during the repetition process, the process proceeds to step S10007.

[0083] (vi) Step S10006 The rejection determination unit 202 returns "rejection" as an output, and the process ends.

[0084] (vii) S10007 The rejection determination unit 202 returns "accepted" as an output, and the process ends.

[0085] 11 is a diagram showing a concept of the rejection determination process performed in steps S10002 to S10006. As shown in FIG. 11, in the algorithm example 1, the parameter value q to be determined is a parameter value q other than the target parameter value q ｕｎ，ｉ The closer the parameter value q to be judged is, the higher the probability of rejecting the parameter value q. This probability is expressed by the probability density function f ｉ More specifically, the parameter value q is calculated using the unintended parameter value q ｕｎ，ｉ If the parameter value q matches the parameter value q, the parameter value q is rejected. ｕｎ，ｉ The further away from , the smaller the probability of rejection. ｕｎ，ｉ The further away from X ｉ The probability that q falls within the range 1101 in FIG. 11 decreases (q becomes less likely to be rejected).

[0086] (2) Algorithm Example 2: Rejection Decision Processing Algorithm 1 Using Euclidean Distance FIG. 12 is a flowchart for explaining algorithm example 2 of the rejection decision function.

[0087] (i) Step S12001 The rejection determination unit 202 of the parameter value generation unit 200 receives as input a parameter value q to be determined and a list Q of parameter values ​​that are not of interest. ｕｎThen, the rejection determination unit 202 receives the list of unstable parameter values ​​Q ｕｎ Take one parameter value from q ｕｎ，ｉ Then, the following processing from step S12002 to step S12004 is performed.

[0088] (ii) Step S12002 The rejection determination unit 202 determines whether the parameter value q to be determined and the parameter value q outside the target range are correct. ｕｎ，ｉ The Euclidean distance d ｉ In addition to the Euclidean distance, the standard Euclidean distance, the Mahalanobis distance, the Manhattan distance, the Chebyshev distance, or the Minkowski distance can be used.

[0089] (iii) Step S12003 The rejection determination unit 202 generates a random number X that follows a continuous uniform distribution U(0,1). ｉ Generates.

[0090] (iv) Step S12004 The rejection determination unit 202 ｉ / d ｍａｘ <X ｉ Determine whether or not d ｉ / d ｍａｘ <X ｉ If it is determined that the above is true (Yes in step S12004), the process proceeds to step S12005. ｉ / d ｍａｘ ≧X ｉ (No in step S12004), the rejection determination unit 202 ｕｎ to another parameter value q ｕｎ，ｉ＋１ Then, the process returns to step S12002. ｍａｘ The specific value of is assumed to be part of the rejection decision function and is not received as an external input. ｍａｘ The kinetic model creation and metabolic simulation system 40 and the rejection judgment function may be configured and defined so that the value of can be included in the input information 001 and input to the rejection judgment function.

[0091] As described above, by the processing from step S12002 to step S12004, q ｕｎ，ｉThis series of processes is performed using Q ｕｎ If the process proceeds to step S12005 due to the conditional branch in step S12004 during the repetition, the parameter values ​​that have not yet been used in the process are ｕｎ If the process does not proceed to S12005 even once due to the conditional branch in step S12004 during the repetition process, the process proceeds to step S12006.

[0092] (v) Step S12005 The rejection determination unit 202 returns "rejection" as an output, and the process ends.

[0093] (vi) Step S12006: The rejection determination unit 202 returns "accepted" as an output, and the process ends. FIG. 13 is a diagram showing a concept of the rejection determination process performed by the processes of steps S12002 to S12005. As shown in FIG. 11, in the algorithm example 2, when the parameter value q to be determined is an unintended parameter value q ｕｎ，ｉ The closer the parameter value q to be judged is, the higher the probability of rejecting the parameter value q. ｕｎ，ｉ Distance d ｉ That is, the probability that the parameter value q to be judged is rejected is calculated using the following formula: ｕｎ，ｉ 11, more specifically, when the parameter value q to be determined is an unintended parameter value q ｕｎ，ｉ If it matches, d ｉ / d ｍａｘ = 0, so X ｉ When takes a random number between 0 and 1, X ｉ The parameter q is rejected regardless of the value of . ｉ / d ｍａｘ If is close to 1, X ｉ If q does not take a value close to 1, the parameter q is accepted (not rejected) with a high probability.

[0094] (3) Algorithm Example 3: Rejection Decision Processing Algorithm 2 Using Euclidean Distance FIG. 14 is a flowchart for explaining algorithm example 3 of the rejection decision function.

[0095] (i) Step S14001 The rejection determination unit 202 of the parameter value generation unit 200 receives as input a parameter value q to be determined and a list Q of parameter values ​​that are not of interest. ｕｎ (List of unstable parameter values ​​306 in the first embodiment). Then, the rejection determination unit 202 receives Q ｕｎ Take one parameter value from q ｕｎ，ｉ Then, the following steps S14002 and S14003 are performed.

[0096] (ii) Step S14002 The rejection determination unit 202 determines whether q and q ｕｎ，ｉ The Euclidean distance d ｉ In addition to the Euclidean distance, the standard Euclidean distance, the Mahalanobis distance, the Manhattan distance, the Chebyshev distance, or the Minkowski distance can be used.

[0097] (iii) Step S14003 The rejection determination unit 202 ｉ <d ｔｈｒｅｓｈ Determine whether or not d ｉ <d ｔｈｒｅｓｈ If it is determined that the value is d (Yes in step S14003), the process proceeds to step S14004. ｉ ≧d ｔｈｒｅｓｈ (No in step S14003), the rejection determination unit 202 ｕｎ to another parameter value q ｕｎ，ｉ＋１ Then, the process returns to step S14002. ｔｈｒｅｓｈ The specific value of is assumed to be configured as part of the rejection decision function and not received as an external input. ｔｈｒｅｓｈ The kinetic model creation and metabolic simulation system 40 and the rejection judgment function may be configured and defined so that the value of can be included in the input information 001 and input to the rejection judgment function.

[0098] As described above, by the processes in steps S14002 and S14003, q ｕｎ，ｉ This series of processes is performed using Q ｕｎ If the process proceeds to S14004 due to the conditional branch in S14003 during the repetition, the parameter value that has not yet been used in the process is Q ｕｎ If the process does not proceed to step S14004 even once due to the conditional branch in step S14003 during the repetition process, the process proceeds to step S14005.

[0099] (iv) S14004 The rejection determination unit 202 returns "rejection" as an output, and the process ends.

[0100] (v) S14005 The rejection determination unit 202 returns "accepted" as an output, and the process ends.

[0101] 15 is a diagram showing a concept of the rejection determination process performed in steps S14002 to S14004. As shown in FIG. 15, in the algorithm example 3, the parameter value q to be determined is the parameter value q ｕｎ，ｉ Distance d from ｔｈｒｅｓｈ If the parameter value q is within the range of q, the parameter value q to be judged is rejected.

[0102] Next, an example will be described in which a set of stable parameter values ​​was obtained by steps S6001 to S6008 (see FIG. 6 ) according to the embodiment described above. Note that the kinetic model having unknown parameters included in input information 001 was obtained by modifying the kinetic model of central carbon metabolism of Escherichia coli (https: / / doi.org / 10.1016 / j.jbiosc.2017.09.005), whose parameters had already been determined, with the following modifications (i) to (iii):

[0103] (i) The dependent variables, external glucose concentration, external acetate concentration, and protein concentration, were fixed at steady-state values ​​in continuous culture at a dilution rate of 0.2 and treated as constants. (ii) Dependent variables, rate equations, and parameters that became unnecessary in (i) were deleted. (iii) The k of phosphoenolpyruvate synthase was ｃａｔ is a parameter representing Ｐｐｓ＿ｃａｔ was treated as an unknown parameter. Sampling and rejection were performed using k Ｐｐｓ＿ｃａｔ The common logarithm was used as the target, and when setting it in the kinetic model, it was converted from the common logarithm to an antilogarithm.

[0104] The target set size N included in the input information 001 was set to 100. The proposal distribution included in the input information 001 was defined as a normal distribution with a mean of 8.06 and a variance of 2.55 (standard deviation σ=1.60).

[0105] The rejection decision function included in the input information 001 is defined as in the above algorithm example 1. However, the variance σ ２ is 2.55 x 10 －２ In addition, the threshold value for determining whether or not the steady state is present in step S6006 is set to 10 －６ It was decided.

[0106] 16 to 18 are diagrams showing the results of collecting stable parameter values ​​in this embodiment and, for comparison, in a conventional EM that does not use the empirical rejection method. In Figures 16 to 18, "with empirical rejection method" refers to the results of this embodiment, and "without empirical rejection method" refers to the results when the rejection determination function is not used in step S6003 of this embodiment and a "accept" decision is always made, i.e., when conventional EM is performed.

[0107] 16 is a graph showing how stable parameter values ​​are collected over the course of calculation time in this example. As shown in FIG. 16, it was confirmed that the introduction of the empirical rejection method reduced the calculation time required to collect 100 stable parameter values ​​to two-thirds or less.

[0108] FIG. 17 is a graph showing how unstable parameter values ​​are recorded during the stable parameter value collection process of FIG. 16 according to this embodiment. As shown in FIG. 17 , it was confirmed that the introduction of empirical rejection reduced the number of unstable parameter values ​​recorded until 100 stable parameter values ​​were collected, from 119 to 26, a reduction of more than one-quarter. This indicates an improvement in the probability that parameter values ​​used for stability assessment are stable. Furthermore, the reason for the reduction in calculation time using empirical rejection, as shown in FIG. 17 , is believed to be because stability assessments were not performed the number of times equal to the difference in the number of unstable parameter values ​​recorded. In other words, unstable parameter values ​​accumulate over time, and parameter values ​​close to these unstable parameter values ​​are discarded before being recorded, thereby reducing calculation time.

[0109] 18 is a graph showing the final distribution of parameter values ​​determined to be stable and unstable in FIGS. 16 and 17 according to this example. As shown in FIG. 18, 97 unstable parameter values ​​(10 ７ ~10 １０ When the marker was introduced, the number of markers decreased to 15.

[0110] On the other hand, looking at the distribution of stable parameter values ​​(distribution on the left side of Figure 18), it can be seen that by introducing the empirical rejection method, the range of possible stable parameter values ​​is narrowed, and the distribution of parameter values ​​approaches the original values ​​without falling into a local optimum.

[0111] <Summary> In this embodiment, in EM where the parameter value processing step 12 involves numerical integration, the total calculation time is reduced by preventing the collected target parameter values ​​from falling into a local solution and increasing the probability that the target parameter values ​​are provided to the parameter value processing step 12. EM is a known technique, but this embodiment is characterized in that it introduces a new technique called empirical rejection.

[0112] The functions of the present embodiment can also be realized by software program code. In this case, a storage medium on which the program code is recorded is provided to a system or device, and the computer (or CPU or MPU) of the system or device reads the program code stored in the storage medium. In this case, the program code read from the storage medium itself realizes the functions of the above-described embodiment, and the program code itself and the storage medium on which it is stored constitute the present disclosure. Examples of storage media for providing such program code include flexible disks, CD-ROMs, DVD-ROMs, hard disks, optical disks, magneto-optical disks, CD-Rs, magnetic tapes, non-volatile memory cards, and ROMs.

[0113] Furthermore, an operating system (OS) running on a computer may perform some or all of the actual processing based on instructions in the program code, and the functions of the above-described embodiments may be realized by this processing.Furthermore, after the program code is read from a storage medium and written to memory on the computer, a CPU of the computer may perform some or all of the actual processing based on instructions in the program code, and the functions of the above-described embodiments may be realized by this processing.

[0114] Furthermore, the program code of the software that realizes the functions of the embodiments and examples may be distributed via a network and stored in a storage means such as a hard disk or memory of the system or device, or in a storage medium such as a CD-RW or CD-R, so that when in use, the computer (or CPU or MPU) of the system or device reads and executes the program code stored in the storage means or storage medium.

[0115] The processes and techniques described herein are not inherently related to any specific device and can be implemented by a combination of components. Various types of general-purpose devices can also be added. A dedicated device may be constructed to perform the functions of the present embodiment and examples. Various functions can also be formed by appropriately combining multiple components disclosed in the present embodiment and examples. For example, some components may be omitted from all the components shown in the embodiments and examples, or components from different examples may be appropriately combined.

[0116] Although specific examples are described in this disclosure, they are in all respects for the purpose of explanation (understanding the technology of the present disclosure) and not for the purpose of limitation. Those skilled in the art will recognize that there are many combinations of hardware, software, and firmware suitable for implementing the technology of the present disclosure. For example, the software described can be implemented in a wide variety of programming or scripting languages, such as assembler, C / C++, Perl, Shell, PHP, Java (registered trademark), etc.

[0117] Furthermore, in the above-described embodiment, the control lines and information lines are those that are considered necessary for the explanation, and not all control lines and information lines in the product are necessarily shown. All components may be interconnected.

[0118] In addition, other implementations of the present disclosure will be apparent to those skilled in the art from consideration of the present embodiments and examples. The specification and examples are exemplary only, with the scope and spirit of the present disclosure being indicated by the following claims.

[0119] 11...sampling process, 12...parameter value processing process, 13...termination determination process, 21...training data, 22...training process, 23...generation task execution process, 31...rejection determination process, 32...unintended parameter value, 40...kinetic model creation and metabolic simulation system, 41...control device, 42...storage device, 43...input device, 44...output device, 001...input information, 100...input receiving unit, 200...parameter value generation unit, 201...sampling unit, 202...rejection determination unit, 300...parameter value processing unit, 301...parameter setting unit, 302...Numerical integration unit, 303...Steady state determination unit, 304...Linear stability analysis unit, 305...Stability recording unit, 306...List of unstable parameter values, 400...Termination control unit, 401...Termination determination unit, 402...Parameter value output unit, 403...List of stable parameter values, 500...Perturbation effect prediction unit, 501...Parameter setting unit, 502...MCA unit, 503...Statistical processing unit, 504...Prediction result output unit, 505...Perturbation effect prediction result, 700...Input screen, 800...Stable parameter value list output screen, 900...Perturbation effect prediction result output screen

Claims

1. A kinetic model generation system that collects target parameter values ​​for a kinetic model having unknown parameters, and generates the kinetic model by setting the collected target parameter values ​​as values ​​of the unknown parameters, comprising: a storage device that stores a program for generating a set of target parameter values ​​representing values ​​of the unknown parameters that contribute to the kinetic model having desired properties, based on untarget parameter values ​​representing values ​​of the unknown parameters that do not contribute to the kinetic model having desired properties; and a processor that reads and executes the program from the storage device, wherein the processor performs: a sampling process that generates candidate parameter values ​​to be set for the unknown parameters by generating random numbers from a given proposal distribution; a rejection determination process that determines whether or not to accept the candidate parameter values ​​based on the untarget parameter values; and a parameter value process that determines whether or not the kinetic model has the desired properties when the candidate parameter values ​​accepted by the rejection determination process are applied to the kinetic model.

2. A kinetic model generation system according to claim 1, wherein the processor, in the parameter value processing, performs numerical integration on the kinetic model to which the candidate parameter value accepted by the rejection determination processing is applied, thereby determining whether the kinetic model to which the candidate parameter value is applied has the desired properties.

3. A kinetic model generation system as set forth in claim 2, wherein the processor, in the parameter value processing, determines that the kinetic model is in a steady state and has the desired properties when the absolute value of the relative rate of change of the metabolite concentration in the state calculated based on the state of the kinetic model in the final step obtained by the numerical integration and the kinetic model is smaller than a predetermined threshold value.

4. A kinetic model generation system according to claim 3, wherein the processor, in the parameter value processing, further performs a linear stability analysis on the kinetic model in the steady state, determines whether the steady state is asymptotically stable, and determines that the asymptotically stable kinetic model has the desired properties.

5. A kinetic model generation system according to claim 1, wherein the processor further executes a process of adding the candidate parameter value to the unintended parameter value when the parameter value process determines that the kinetic model does not have the desired properties.

6. A kinetic model generation system as claimed in claim 1, wherein the processor, in the rejection determination process, calculates the value of a probability density function when the value of a random variable in a probability distribution that includes the undesired parameter value as a parameter is the candidate parameter value generated in the sampling process, and determines whether or not to accept the candidate parameter value based on the value of the probability density function.

7. A kinetic model generation system according to claim 6, wherein the processor accepts parameters other than the untarget parameter value input from outside, and calculates the value of the probability density function according to the probability distribution of the accepted parameters.

8. A kinetic model generation system according to claim 6, wherein the probability distribution is any one of a normal distribution, a log-normal distribution, a uniform distribution, a gamma distribution, a triangular distribution, and a mixed distribution including at least two of them.

9. A kinetic model generation system according to claim 6, wherein the processor determines whether or not to accept the candidate parameter values ​​generated in the sampling process based on the magnitude relationship between the value of the probability density function and a separately generated uniform random number.

10. A kinetic model generation system according to claim 1, wherein the processor, in the rejection determination process, calculates the distance between the candidate parameter value generated in the sampling process and the unintended parameter value, and determines whether or not to accept the candidate parameter value based on the distance.

11. A kinetic model generation system according to claim 10, wherein the processor calculates the distance as any one of Euclidean distance, standard Euclidean distance, Mahalanobis distance, Manhattan distance, Chebyshev distance, and Minkowski distance.

12. A kinetic model generation system according to claim 11, wherein the processor determines whether or not to accept the candidate parameter values ​​generated in the sampling process based on the magnitude relationship between the distance and a predetermined threshold value.

13. A kinetic model generation system according to claim 12, wherein the processor receives the threshold value input from an external source and determines whether the received threshold value is larger or smaller than the distance.

14. A kinetic model generation system according to claim 1, wherein the processor further executes an end determination process that determines whether to end a series of processes including the sampling process, the rejection determination process, and the parameter value process, or to start again from the sampling process.

15. A kinetic model generation system according to claim 14, wherein the processor terminates the series of processes based on the total number of target parameter values ​​reaching a predetermined value in the termination determination process.

16. A kinetic model generation system according to claim 1, wherein the processor accepts the desired properties of the kinetic model or the proposed distribution input from outside, and executes the program.

17. The kinetic model generation system according to claim 1, wherein the desired property of the kinetic model is stability.

18. A kinetic model generation system according to claim 1, wherein the processor determines whether the kinetic model has the desired properties based on the value of an evaluation function calculated using actual data and states or state transitions obtained by performing numerical integration on the kinetic model.

19. The kinetic model generation system according to claim 18, wherein the evaluation function is one of mean square error, root mean square error, cross entropy error, mean absolute error, or squared logarithmic error.

20. A kinetic model generation system according to claim 1, wherein the processor determines whether the kinetic model has the desired properties based on the time at which the kinetic model reaches a steady state in the results obtained by performing numerical integration on the kinetic model.

21. A kinetic model generation system according to claim 1, wherein the processor executes a process of outputting at least one of the set of target parameter values ​​and the set of kinetic models obtained by applying parameter values ​​included in the set of target parameter values.

22. A metabolic control analysis system comprising the kinetic model generation system of claim 1, wherein the processor performs metabolic control analysis using the set of target parameter values ​​in the kinetic model.

23. A metabolic control analysis system according to claim 22, wherein the processor executes a process for outputting the results of the metabolic control analysis in a viewable and / or usable format.

24. A metabolic control analysis system according to claim 22, wherein the processor statistically processes the results of the metabolic control analysis.

25. A metabolic control analysis system according to claim 24, wherein the processor executes a process for outputting the results of the statistical processing in a viewable and / or usable format.

26. A kinetic model generation method for collecting target parameter values ​​for a kinetic model having unknown parameters, and generating the kinetic model by setting the collected target parameter values ​​as values ​​of the unknown parameters, the method comprising: a processor generating a set of target parameter values ​​representing values ​​of the unknown parameters that contribute to the kinetic model having desired properties, based on non-target parameter values ​​representing values ​​of the unknown parameters that do not contribute to the kinetic model having desired properties; wherein generating the set of target parameter values ​​comprises: the processor executing a sampling process to generate candidate parameter values ​​to be set for the unknown parameters by generating random numbers from a given proposal distribution; the processor executing a rejection decision process to determine whether or not to accept the candidate parameter values, based on the non-target parameter values; and the processor executing a parameter value process to determine whether or not the kinetic model has the desired properties when the candidate parameter values ​​accepted by the rejection decision process are applied to the kinetic model.

27. A program for causing a computer to execute the kinetic model generation method according to claim 26.

Images