Method for realizing efficient solution processing for chemical principal equation based on variational auto-encoder and quantization tensor sequence

Abstract

The invention relates to a method for realizing efficient solution processing for a chemical main equation based on a variational auto-encoder and a quantitative tensor sequence, which comprises the following steps of: S1, carrying out random simulation on a specific biochemical reaction by utilizing a random simulation algorithm to obtain the number of target molecules, and calculating the probability distribution of the target molecules as a data set; s2, inputting the initial probability distribution of the target molecule into a trained variational auto-encoder, mapping the probability distribution to an equivalent tendency function by using the variational auto-encoder, and performing iteration to obtain an equivalent reaction tendency function of an intermediate reaction; and S3, inputting different biochemical reaction kinetic parameters into the trained variational auto-encoder to obtain a probability distribution solution of the chemical main equation. By adopting the method for efficiently solving the chemical main equation based on the variational auto-encoder and the quantitative tensor sequence, the complexity of solving the chemical main equation of a biochemical reaction system is reduced, the sample size of training is reduced, the training time is shortened, and a powerful tool is provided for analyzing a random biochemical system.

Description

Technical Field

[0001] The present invention relates to the technical field of stochastic kinetic modeling of biochemical reactions, in particular to the efficient solution of the chemical master equation, and specifically refers to a method for efficiently solving and processing the chemical master equation based on a variational autoencoder and a quantized tensor sequence. Background Art

[0002] Studying biochemical reaction kinetics at the mesoscopic scale is crucial for understanding complex processes in chemical engineering and biological systems. Reactions at the mesoscopic scale include crystallization reactions, polymerization reactions, and intracellular biochemical reactions. Research at this scale is the key to revealing the dynamics of chemical engineering and biological systems. The mesoscopic scale lies between the macroscopic and microscopic scales. At this scale, reactants are regarded as molecular ensembles, focusing on the number of molecules while ignoring the internal structure of the molecules. By simplifying each molecule into a sphere, the mesoscopic scale method can capture important phenomena ignored by macroscopic models and can also perform long-term dynamic predictions that are difficult to achieve with microscopic models. Although this method ignores microscopic structural details, in many reaction systems, due to the small number of molecules and the inherent random interactions between molecules, the fluctuations in the number of molecules become significant. Therefore, probability distribution becomes the core tool for describing mesoscopic reaction kinetics.

[0003] The Chemical Master Equation (CME) is exactly the mathematical tool used to describe this randomness and discreteness. The CME consists of a set of differential-difference equations and can effectively capture the randomness and discreteness in the reaction process. The CME is particularly suitable for reaction systems at the mesoscopic scale, where the number of molecules is small and continuous concentration approximation cannot be used. By considering the randomness and discreteness of chemical reactions, the CME can characterize the fluctuation characteristics of the number of molecules. In the CME, the number of molecules is discrete, which is different from traditional continuous concentration models. Since chemical reactions are inherently random, the CME captures the probability evolution process of each possible state.

[0004] To date, many methods have been proposed to solve the CME, among which the Finite State Projection (FSP) method and the Gillespie algorithm (SSA) are particularly prominent. However, due to the curse of dimensionality, these methods have encountered significant challenges in dealing with large and complex systems. For example, to fully simulate a gene expression system with M species (each with a maximum number of N), it is usually necessary to solve (N + 1)^M ordinary differential equations. This means that when solving the CME, the probability distributions of these state combinations and their changes over time must be considered simultaneously. As the dimensionality of the state space increases, the number of equations in the system grows exponentially, greatly increasing the computational difficulty. To address this challenge, researchers have proposed various approximation methods, such as continuous approximation, time-scale separation, and linear noise approximation. However, these methods usually rely on a series of restrictive assumptions that may not be applicable to a wide range of reaction systems.

[0005] In the past decade, the rapid development of neural networks has promoted technological progress in fields such as deep learning, machine learning, and data science. These technologies have become important tools in modern scientific research and engineering practice, and their influence has almost covered all fields, including chemical engineering, physics, materials science, biology, etc. In the field of chemical engineering, one of the core roles of neural networks is to utilize their powerful modeling and prediction capabilities to solve the curse of dimensionality problem in complex kinetic systems.

[0006] As a tensor decomposition technique, the Quantized Tensor Train (QTT) method has received extensive attention in recent years, especially in the fields of machine learning and data science. Traditional high-dimensional tensor storage requires a large amount of memory and disk space. Through QTT decomposition, tensors can be stored as multiple low-rank matrices, and further reduction of the required storage space can be achieved through quantization. For large-scale data sets, QTT can significantly reduce the storage complexity. The quantized QTT not only reduces the storage requirements but also reduces the computational amount. Especially in matrix operations and tensor operations, the demand for computing resources is greatly reduced. This makes QTT particularly suitable for application scenarios with limited computing resources. The QTT method can handle high-dimensional data and avoid the "curse of dimensionality" by decomposing high-dimensional tensors into the product of multiple low-rank matrices. This is very effective for many practical problems, such as high-dimensional data compression, dimensionality reduction, and pattern recognition.

[0007] In summary, neural networks and quantized tensor trains are of great significance and broad prospects in solving high-dimensional chemical master equations. Through continuous optimization and innovation, related technologies can provide more reliable and efficient solutions for the efficient solution of chemical master equations. Summary of the Invention

[0008] The object of the present invention is to overcome the above-mentioned disadvantages of the prior art and provide a method for efficiently solving and processing chemical master equations based on variational autoencoders and quantized tensor sequences. Although the application prospects of neural networks in chemical engineering are broad, there are still some challenges. First, the "black box" nature of neural networks limits their scientific interpretability; second, neural networks require a large amount of labeled data for training, and in biochemical reaction systems, the cost of obtaining experimental data is often high; the lack of physical interpretability of the model, the strong dependence on large-scale high-quality data, and the high computational resource requirements during training are all problems that need to be solved currently. Training a neural network-based model is both time-consuming and computationally intensive. Therefore, a key challenge is to develop a neural network that can be generalized to different systems. However, the prediction accuracy of a model trained on a specific reaction system is not good when applied to other systems. A feasible method is to derive equivalent reactions from a set of intermediate reactions, which can reduce the dimensionality and improve the computational efficiency. However, there are still few methods for constructing such effective reactions in biochemical systems.

[0009] To solve this problem, the present invention uses a variational autoencoder (VAE) as a powerful tool for identifying effective reactions in complex systems, thereby significantly simplifying the CME. The VAE effectively learns the latent representation by encoding the input data into a low-dimensional space and capturing the key features of the data. Its encoding-decoding process allows capturing complex patterns and dependencies in the data, making it very suitable for identifying equivalent reactions. At the same time, by utilizing the data compression ability of the quantization tensor sequence, the complex CME is simplified again, thus solving the curse of dimensionality. The method for efficiently solving the chemical master equation based on the variational autoencoder and the quantization tensor sequence includes: using the stochastic simulation algorithm to perform stochastic simulation on specific biochemical reactions to obtain the number of target molecules and calculate the probability distribution of the target molecules as a data set; inputting the initial probability distribution of the molecules into the trained variational autoencoder, iteratively obtaining the equivalent reaction propensity function of the intermediate reaction, and then transforming the chemical master equation into the form of a quantization tensor sequence for solution. The training steps of the variational autoencoder include: dividing the entire reaction set into two non-overlapping subsets, respectively representing the reactions of interest or not; according to the division result of the entire reaction system, determining the equivalent propensity function that the variational autoencoder needs to learn and obtaining a simplified chemical master equation; inputting the initial probability distribution of the molecules of interest into the encoder, dividing the output into two parts, which are used as the mean and variance of a Gaussian distribution; performing a reparameterization operation on the output of the encoder to obtain the input of the decoder, making it possible to calculate the gradient; inputting the result of the parameterization into the decoder, adding a bias term to the obtained output, and finally obtaining the output of the variational autoencoder, that is, the propensity function of the equivalent reaction; using the output of the variational autoencoder to calculate the probability distribution at the next time point, repeating the above steps until the termination time of the reaction is reached; after obtaining the equivalent reaction propensity function of the intermediate reaction, obtaining a simplified chemical master equation and converting the simplified chemical master equation into the form of a quantization tensor sequence, thereby efficiently solving the complex chemical master equation; calculating the corresponding loss function by solving the chemical master equation to obtain the result, and repeating the above steps for training. Through the method of the variational autoencoder and the quantization tensor sequence, the present invention realizes the improvement of the solving speed of the chemical master equation and can be applied to different chemical reaction systems, having certain practical value.

[0010] Preferably, in some examples of the present invention, the step of generating the data set includes: using the Stochastic Simulation Algorithm (SSA) to obtain the data set required for training. SSA is a numerical method for simulating discrete-time stochastic processes, especially suitable for modeling chemical reaction networks, and is widely used in bioinformatics, systems biology, and biochemical kinetics. Its core is to simulate the occurrence of chemical reactions through a discretized stochastic process, rather than relying on differential equations like traditional continuous kinetic models. The advantage of this algorithm is that it can accurately capture the reaction fluctuations in the case of low molecule numbers, especially important when dealing with small-scale molecular systems. For different biochemical reactions, various improved stochastic simulation algorithms have been developed to better simulate these special biochemical reactions, such as the Delay Rejection Algorithm, the Delay Direct Method Algorithm, the Delay Modified Next Reaction Method Algorithm, and so on. This algorithm is an accurate stochastic simulation method that can fully restore the stochastic behavior of the system when the number of molecules is small, especially suitable for dealing with small-scale reaction systems with noise. It is applicable to various complex chemical reaction networks, including models with arbitrary complex reactions and any number of substances, and can simulate the time evolution process of the system without assuming that the system reaches a steady state or equilibrium.

[0011] Preferably, the initial probability distribution of the molecules is input into the trained variational autoencoder, and the equivalent reaction propensity function of the intermediate reaction is obtained iteratively. Then, the chemical master equation is transformed into a tensor sequence form to efficiently solve the chemical master equation. The training steps of the variational autoencoder include:

[0012] Step S21, considering a reaction system containing molecules X i (i = 1,…, P) and R reactions, classifying the indices of these molecules into a set Dividing the entire reaction set into two non-overlapping subsets A and B, with R A and R B reactions respectively. The reaction kinetics in set A is the main concern, while the reaction kinetics in set B is the secondary concern and contains many unimportant intermediate reactions. Divide all molecules into three categories: Set contains the indices of molecules involved in both set A and B, contains the indices of molecules only related to set A, contains the indices of molecules only related to set B;

[0013] Step S22: According to the partitioning result of the entire reaction system, determine the equivalent propensity function that the variational autoencoder needs to learn, and obtain a simplified chemical master equation;

[0014] Step S23: Input the initial probability distribution of the molecule of interest into the encoder, and divide the output into two parts, which are used as the mean and variance of a Gaussian distribution;

[0015] Step S24: Perform a reparameterization operation on the output of the encoder to obtain the input of the decoder, making it possible to calculate the gradient;

[0016] Step S25: Input the parameterized result into the decoder, add a bias term to the obtained output, and finally obtain the output of the variational autoencoder, which is the equivalent reaction propensity function;

[0017] Step S26: Use the output of the variational autoencoder to calculate the probability distribution at the next time point, and repeat the above steps until the termination time of the reaction is reached;

[0018] Step S27: After obtaining the equivalent reaction propensity function of the intermediate reaction, obtain a simplified chemical master equation, and convert the simplified chemical master equation into a form of a quantization tensor sequence, so as to efficiently solve the complex chemical master equation;

[0019] Step S28: The loss function consists of the reconstruction error and the KL divergence. Input the result obtained by solving the chemical master equation back into the encoder, and the obtained output is used to calculate the KL divergence. Calculate the gap between the result obtained by solving the chemical master equation and the probability distribution obtained by the stochastic simulation algorithm as the reconstruction error. Repeat steps S22 - 28 until the overall loss function of the network reaches the preset target, and finally obtain the trained variational autoencoder model;

[0020] Preferably, by inputting different biochemical reaction kinetic parameters into the trained variational autoencoder, the chemical master equation of this biochemical reaction can be efficiently solved to obtain the probability distribution solution of the chemical master equation.

[0021] Preferably, in some embodiments of the present invention, step S21 includes:

[0022] Partition the entire reaction set into two non - overlapping subsets A and B, which contain R A and R B reactions respectively. The reaction kinetics in set A is the main concern, while the reaction kinetics in set B is the secondary concern. In many cases, some molecules in set A also participate in the reactions in set B, and some molecules involved in set A are the products of the reactions in set B. Divide the indices of all molecules (P kinds of molecules) into three categories: H is the number of molecular species related to sets A and B, M - H + 1 is the number of molecular species related only to set A, and P - M + 1 is the number of molecular species related only to set B. M, H, and P need to be set manually according to different reaction systems. Then the set contains the indices of the molecules related to sets A and B, contains the molecules related only to set A, contains the molecules related only to set B. Then the entire reaction system can be expressed as:

[0023]

[0024] In this reaction system, s ir and s ′ ir are both non - negative integers, representing the quantities of reactants and products respectively. The reaction coefficient matrix S is defined as:

[0025] S ir = s ′ ir - s ir (2)

[0026] To solve the curse of dimensionality and reduce the computational complexity, a single equivalent reaction is used to replace the reactions in set B. Since this method ignores the intermediate reactions in set B in reaction (1), only the reactants and products of the first reaction (r = R A + 1) and the last reaction (r = R) in set B need to be concerned about. This reaction can accurately capture their influence on the reaction kinetics in set A. Such a reaction can be expressed as:

[0027]

[0028] For example, assume that set B contains the following three reactions X1 → X2 + X3, X2 + X3 → X4, and assume that the kinetics of X1 needs to be solved. Then the specific form of the equivalent reaction (3) is

[0029] Preferably, in some embodiments of the present invention, step S22 includes:

[0030] Assume that the equivalent propensity function of the equivalent reaction (3) is f eff , and the reduced chemical master equation corresponding to the reaction system jointly composed of the reactions in set A and the equivalent reaction (3) can be given by the following formula:

[0031]

[0032] The original reaction (1) contains P types of molecules. After simplification by this method, the number of molecule types is reduced to M, where represents the number of molecules of M types, where is the probability that the system is in state at time t, where n i is the number of molecules of molecule X i . S r is the r-th column of matrix S. f r (n) is the propensity function for the r-th reaction to occur when the system is in state n, and can be expressed as

[0033]

[0034] where k r is the reaction rate of the r-th reaction. For convenience, Ω = 1 is set.

[0035] Utilizing the universal approximation ability of neural networks, VAE (NN θ ) is used to approximate f eff , where θ are the neural network parameters of VAE, and is represented by the following equation:

[0036]

[0037] Assume N represents the maximum number of all M types of molecules. Therefore, equation (4) can be concisely expressed as:

[0038]

[0039] where P(t) is a vector containing all possible states, ranging from to where for a single type of molecule, P(t) = [P(0,t), …, P(N,t)] T , representing the probability distribution of this species. The state transition matrix A θ (t) is defined as A θ (t) = D + N θ (t), where D contains the propensity function f r in the reactions in set A, while N θ (t) contains the effective propensity NN θ (n,t) related to the reactions in set B.

[0040] When comparing the chemical master equation to be solved for the unsimplified reaction system (1) (containing (N + 1) M equations) with the chemical master equation to be solved for the equivalent reaction system after simplifying the reactions in set B ((N + 1)P The computational advantage becomes evident when calculating the number of ODEs required in the equation (where P >> M). The significant reduction in the number of equations highlights the efficiency of this method. The VAE does not explicitly identify which reactions are intermediate reactions, nor does it determine which reactions should be classified as valid reactions. It automatically calculates the equivalent propensity function, and it is necessary to manually select the molecules of primary interest to indirectly achieve the classification of intermediate reactions.

[0041] Preferably, in some embodiments of the present invention, the step S23 includes:

[0042] Design a suitable encoder E(·; φ), which is a neural network, where · represents the input of the neural network and φ is the parameter of the neural network. The main task of the encoder is to map the input data to a distribution in the latent space. This process can be regarded as an approximate inference method, aiming to infer the probability distribution of the input data in the latent space through the learning of the neural network. The encoder does not directly output a fixed latent variable, but estimates the distribution of the latent variable through the network. The encoder is used to estimate the variational distribution, taking P(t) as the input and generating the variational distribution q(z|P(t); φ). Generally, it is assumed that the variational distribution q(z|P(t); φ) follows a Gaussian distribution with a diagonal covariance matrix:

[0043]

[0044] where I is the identity matrix, μ I and σ I are the mean and standard deviation of the Gaussian distribution, z is the latent variable, and these parameters can be predicted by the encoder:

[0045]

[0046] Preferably, in some embodiments of the present invention, the step S24 includes:

[0047] Since the random variable z is sampled from the variational distribution q(z|P(t); φ), z and φ are not deterministically related. To solve this problem and enable the model to perform backpropagation and optimize the parameters, z can be reparameterized as:

[0048] z = μ I + σ I ⊙ ∈(10)

[0049] where and ⊙ represents the Hadamard product. This reparameterization method converts the random relationship between z and φ into a deterministic relationship, so that the derivative of z with respect to φ can be directly calculated. The variational autoencoder can perform effective sampling in the latent space while maintaining the differentiability of the model, enabling it to be trained using standard gradient descent methods.

[0050] Preferably, in some embodiments of the present invention, the step S25 includes:

[0051] Design a suitable decoder D(·; ψ), which is also a neural network, where · represents the input of the neural network and ψ are the parameters of the neural network. Using the reparameterized result z as the data for the decoder, the output of the variational autoencoder can be obtained:

[0052] D(z; ψ) = p(P(t)|z; ψ) (11)

[0053] To accelerate the training, a set of preset biases can be added to encode the prior knowledge about the reaction kinetics in set B. These biases enhance the model by providing suitable initial parameters for the equivalent propensity function. Therefore, the final output of the VAE model is given by the following formula:

[0054]

[0055] Preferably, in some embodiments of the present invention, the step S26 includes:

[0056] Using formula (4), the equivalent propensity function and the input P(t), it is very simple to perform a feed-forward calculation to predict the distribution P(t + Δt) at subsequent times. For example, the Euler method can be used for the solution:

[0057] P(t + Δt) = P(t) + A θ (t)P(t)Δt. (13)

[0058] This method can predict the distribution at any time t ∈ [0, T] given T.

[0059] Preferably, in some embodiments of the present invention, the step S27 includes:

[0060] The matrix A θ (t) is a high-dimensional matrix and directly representing it would result in huge storage requirements. To effectively represent this matrix, the tensor train (TT) format is used to solve this problem. The TT format decomposes a high-dimensional tensor into multiple low-rank matrices. Suppose there is an N-dimensional matrix A θ (t) of size I1 × I2 × … × I N , it can be represented in the tensor train format as:

[0061] A θ (i1, …, i N ) = A1(i1)·A2(i2)·…·A N (i N ) (14)

[0062] where each A k (i k ) is a matrix of rank r k-1 ×r k .

[0063] The Quantized Tensor Train (QTT) format is an extension of the TT format, aiming to further reduce the rank of tensors, especially suitable for systems with large-scale state spaces. In the QTT format, each dimension of the tensor is rearranged so that the size of each mode becomes a power of 2. In this way, through binary representation, the index of each mode is compressed, further reducing the storage requirements. Suppose there is an N-dimensional matrix A θ (t) of size I1×I2×…×I N , first remap the size I k of each mode to the power of 2, 2L k , and then reshape the tensor into a new tensor Y:

[0064] Y(i1,i2,…,i N ) = A H (i1,i2,…,i N ) (15)

[0065] The size of this new tensor Y is 2×2×…×2, and the index of each mode uses binary representation, so that it can be efficiently stored and solved using the QTT format.

[0066] Preferably, in some embodiments of the present invention, the step S28 includes:

[0067] The loss function of the VAE is the sum of the reconstruction error and the difference between the posterior distribution and the prior distribution in the latent space. It can be expressed as:

[0068]

[0069] where β is the weight of the reconstruction error. represents the Kullback-Leibler (KL) divergence, which is used to quantify the difference between the posterior distribution q(z|P(t); φ) and the prior distribution p(z|P(t); ψ). The reconstruction error is the mean square error (MSE) between the obtained probability distribution P(t) obtained by solution and the reference distribution generated by the exact stochastic simulation algorithm (SSA):

[0070]

[0071] where is the set of time points for collecting training data, Nshots is the total number of selected time points. The total KL divergence of all selected time points is given by the following formula

[0072]

[0073] where J is the dimension of the latent space z. Here, μ I = [μ1, …, μ J T and σ I = [σ1, …, σ J T are the mean and standard deviation vectors of the Gaussian distribution.

[0074] By adopting the method for efficiently solving the chemical master equation based on the variational autoencoder and the quantization tensor sequence of the present invention, the complexity of solving the chemical master equation of the biochemical reaction system is reduced, the sample size and training time for training are decreased, and the trained variational autoencoder can predict the kinetics of unseen time points, different reaction structures, and molecules not included in the training data, providing a powerful tool for analyzing stochastic biochemical systems. Brief Description of the Drawings

[0075] Figure 1 Shows a schematic diagram of how a complex biochemical reaction system is divided into two subsets according to the idea of the present invention, and an unimportant intermediate reaction is approximated as an equivalent reaction.

[0076] Figure 2 Shows a schematic diagram of the variational autoencoder network structure and the solution training process in the method for efficiently solving the chemical master equation based on the variational autoencoder and the quantization tensor sequence provided according to some embodiments of the present invention.

[0077] Figure 3 Shows a schematic diagram of the solution results of the chemical master equation applied to three common gene expression models according to the present invention.

[0078] Figure 4 Shows a schematic diagram of the efficiency comparison with other methods on a more complex gene expression model according to the present invention.

[0079] Figure 5 Shows a schematic diagram of the extrapolation method and the extrapolation effect with different reaction kinetic parameters according to the present invention.

[0080] Figure 6 Shows a schematic diagram of the extrapolation method and the extrapolation effect for different reaction structures according to the present invention. Detailed Description of the Embodiments

[0081] ​​In order to more clearly describe the technical content of the present invention, further description is given below in conjunction with specific embodiments.

[0082] Before describing in detail embodiments according to the present invention, it should be noted that, hereinafter, the terms "comprises", "includes" or any other variations are intended to cover non-exclusive inclusion, whereby a process, method, article or apparatus comprising a series of elements includes not only these elements, but also other elements not explicitly listed or inherent to such process, method, article or apparatus.

[0083] See also Figure 1 As shown, the radar-visual fusion target perception method based on the hybrid strategy, wherein the method comprises the following steps:

[0084] First reference Figure 1 , Figure 1 It shows how to divide a complex biochemical reaction system into two subsets and approximate an unimportant intermediate reaction to an equivalent reaction according to the idea of ​​the present invention. Specifically:

[0085] The entire reaction set is divided into two disjoint subsets A and B, each containing R A and R B reactions. In this division, the reaction kinetics in set A are the main focus of the present invention, while the reaction kinetics in set B are of secondary concern. The reactions in set A are usually the core part of the study and may be directly related to the key behaviors of the system; while the reactions in set B have an impact on the dynamics of the overall system, their effect on the reaction kinetics in set A may be small or indirect. In many practical situations, some molecules in set A will also participate in the reactions in set B, and vice versa. For example, some molecules involved in set A may act as reactants in the reactions in set B, or some molecules are products of the reactions in set B, and there is a certain intersection or interaction between the two during the reaction process.

[0086] In order to more clearly represent and manage these molecules and their relationships, all molecules participating in the reaction (a total of P molecules) can be classified according to their roles in the reaction sets A and B. Specifically, the present invention divides the molecular index into three categories: in Represents the complete set of all molecule indices. Contains the indices of molecules that participate in both reactions in set A and reactions in set B; set contains the indices of molecules that participate only in reactions in set A; while set Include the indices of the molecules that only participate in the reactions of set B. Such classification helps to clarify the roles of different molecules in the reaction network, and further provides a clearer framework for the construction and analysis of the model.

[0087] However, since the influence of the reactions in set B on the reaction kinetics of set A is relatively indirect, and the reactions in set B itself may involve a large number of molecules and complex reaction pathways, directly modeling all the reactions in set B may lead to the curse of dimensionality, greatly increasing the computational complexity. Therefore, to solve this problem, the present invention adopts a simplification strategy, that is, by introducing an equivalent reaction to replace all the reactions in set B. This equivalent reaction can accurately capture the influence of the reactions in set B on the reaction kinetics of set A, thereby effectively reducing the computational amount and maintaining the accuracy of the model. In this way, while reducing the computational complexity, it can ensure that the kinetic behavior of the reaction system is accurately described, avoiding problems such as numerical instability or low computational efficiency caused by excessive complication.

[0088] Furthermore, in order to approximate the equivalent propensity function of the equivalent reaction of a series of intermediate reactions, the present invention uses a variational autoencoder to map the probability distribution to the equivalent propensity function, and uses a quantization tensor sequence to obtain the final solution of the chemical master equation. The variational autoencoder network structure and the solution training process are as Figure 2 shown. Specifically:

[0089] First, the initial probability distribution of the molecules of interest is passed as input to the encoder. The encoder transforms these inputs into the parameters of a Gaussian distribution, which include the mean and variance. These parameters are calculated through the feed-forward process of the network, representing the distribution of the molecular states in the latent space. In order to enable gradient update in subsequent calculations, a reparameterization operation is performed on the output of the encoder. Specifically, by appropriately transforming the mean and variance, the network can sample through the standard normal distribution, thereby generating latent variables with differentiable properties. This operation enables the backpropagation algorithm to effectively calculate the gradient during the training process, thereby updating the network parameters.

[0090] Next, the reparameterized latent variables are input into the decoder. The role of the decoder is to map the information in the latent space back to the original data space. Here, the output generated by the decoder is the propensity function of the reaction, and this output represents the reaction rate or reaction probability of the chemical reaction. In order to ensure that the output can stably and accurately reflect the real process of the reaction, a bias term is added to the output of the decoder to adjust the output of the model, making it better fit the actual reaction dynamics. After this adjustment, the final output of the variational autoencoder is the propensity function of the equivalent reaction, reflecting the dynamic evolution of the reaction system.

[0091] Using the output of the variational autoencoder, the probability distribution of the system at the next time point can be calculated. Specifically, based on the state at the current time point and the propensity function of its equivalent reaction, the possible state distribution of the system at the next time point is inferred. This process is continuously repeated, and at each step, the probability distribution at the next moment is calculated through the alternating operations of the encoder and decoder. This process continues until the reaction reaches a predetermined termination time point or until the evolution of the system is completed.

[0092] Through this method, the propensity functions of the equivalent reactions of the intermediate reactions are obtained. Then, these propensity functions of the equivalent reactions are used to simplify the chemical master equation. Traditional chemical master equations usually contain a large number of reaction paths and species, while here, through the propensity functions generated by the variational autoencoder, these complex reaction paths can be simplified into a more compact and efficient form. This simplified chemical master equation can be transformed into a sequence of quantization tensors, thereby being solved more efficiently. This method significantly reduces the computational cost when solving complex reaction systems and improves the solution efficiency at the same time.

[0093] When training the variational autoencoder model, the loss function consists of two main parts: the reconstruction error and the KL divergence. The reconstruction error is used to measure the difference between the model output and the actual reaction data, while the KL divergence is used to measure the difference between the distribution learned in the latent space and the standard normal distribution. During the training process, the results obtained by solving the chemical master equation are re-input into the encoder to generate a new distribution in the latent space, and its KL divergence from the standard normal distribution is calculated. At the same time, the reconstruction error of the model output is calculated, that is, the difference between the probability distribution calculated by the variational autoencoder and the probability distribution obtained by the stochastic simulation algorithm.

[0094] This process continues until the overall loss function of the entire network reaches a preset target, usually meaning that both the reconstruction error and the KL divergence reach a minimum value or are close to the optimal state. Finally, after multiple rounds of iterative training, the variational autoencoder model will be able to accurately simulate the dynamic evolution of chemical reactions and generate effective simplified reaction equations to complete the efficient solution of complex chemical master equations.

[0095] Furthermore, the training process of the entire variational autoencoder network is as follows:

[0096] S1, Using the stochastic simulation algorithm, perform stochastic simulation for a specific biochemical reaction to obtain the number of target molecules, and calculate the probability distribution of the target molecules as the data set;

[0097] S2, Input the initial probability distribution of the molecules into the trained variational autoencoder, and iteratively obtain the propensity functions of the equivalent reactions of the intermediate reactions. The training steps of the variational autoencoder include:

[0098] S21, divide the entire reaction set into two non - overlapping subsets A and B, with R A and R B reactions respectively. The reaction kinetics in set A is the main concern, while the reaction kinetics in set B is a secondary concern, containing many unimportant intermediate reactions. Divide all molecules into three categories: Set contains the indices of molecules involved in both sets A and B, contains the indices of molecules related only to set A, contains the indices of molecules related only to set B;

[0099] S22, according to the partitioning result of the entire reaction system, determine the equivalent propensity function that the variational auto - encoder needs to learn, and obtain the simplified chemical master equation;

[0100] S23, input the initial probability distribution of the molecules of interest into the encoder, and divide the output into two parts, which are used as the mean and variance of a Gaussian distribution;

[0101] S24, perform a re - parameterization operation on the output of the encoder to obtain the input of the decoder, making it possible to calculate the gradient;

[0102] S25, input the parameterized result into the decoder, add a bias to the output, and finally obtain the output of the variational auto - encoder, which is the equivalent reaction propensity function;

[0103] S26, use the output of the variational auto - encoder to calculate the probability distribution at the next time point, and repeat the above steps until the termination time of the reaction is reached;

[0104] S27, after obtaining the equivalent reaction propensity function of the intermediate reactions, obtain the simplified chemical master equation, and convert the simplified chemical master equation into a quantization tensor sequence form to efficiently solve the complex chemical master equation;

[0105] S28, the loss function consists of the reconstruction error and the KL divergence. Re - input the result obtained by solving the chemical master equation into the encoder, and the output is used to calculate the KL divergence. Calculate the gap between the result obtained by solving the chemical master equation and the probability distribution obtained by the stochastic simulation algorithm as the reconstruction error. Repeat steps S22 - S28 until the overall loss function of the network reaches the preset target, and finally obtain the trained variational auto - encoder model;

[0106] S3, input different biochemical reaction kinetic parameters into the trained variational auto - encoder, and then the chemical master equation of this biochemical reaction can be efficiently solved to obtain the probability distribution solution of the chemical master equation.

[0107] Further, the loss function of the network is as follows:

[0108]

[0109] where β is the weight of the reconstruction error. represents the Kullback-Leibler (KL) divergence, which is used to quantify the difference between the posterior distribution q(z|P(t); φ) and the prior distribution p(z|P(t); ψ). The reconstruction error is the mean square error (MSE) between the obtained probability distribution P(t) solved and the reference distribution generated by the exact stochastic simulation algorithm (SSA):

[0110]

[0111] where is the set of time points for collecting training data, and N shots is the total number of selected time points. The total KL divergence of all selected time points is given by the following formula

[0112]

[0113] where J is the dimension of the latent space z. Here, μ I = [μ1, …, μ J T and σ I = [σ1, …, σ J T are the mean and standard deviation vectors of the Gaussian distribution.

[0114] Further, Figure 3 shows the solution results applied to three common gene expression models, the chemical master equation, of the present invention. These three gene expression models are the birth-death model, the burst model, and the telegraph model, respectively.

[0115] The birth-death model describes that genes are transcribed by RNA polymerase II (Pol II), which binds to the promoter and initiates the production of nascent mRNA (N). Nascent mRNA is produced at a constant rate ρ. Pol II molecules move along the gene at a constant speed, which means that after a fixed time interval τ, Pol II will dissociate from the gene. Therefore, nascent mRNA will undergo numerous intermediate reaction steps with the same reaction rate and finally degrade. The present invention regards the reaction set as a delayed reaction. The whole system can be represented as:

[0116]

[0117] ​​The burst model is the same as the birth-death model, except that the binding of Pol II to the promoter occurs in bursts, and the size i follows a geometric distribution b i / (1 + b) i+1 , and the whole system can be expressed as:

[0118]

[0119] where α represents the burst frequency and b represents the average burst size. The burst model is used to explain the phenomenon that gene expression occurs in high-intensity bursts rather than at a steady rate.

[0120] Different from the birth-death model and the burst model, the telegraph model describes gene expression as a process in which a gene can switch between an active (on) and an inactive (off) state. These switches are controlled by two rate constants: the activation rate (σ on ) and the inactivation rate (σ off ). When the gene is in the on state, mRNA is produced at a constant rate, and when the gene is in the off state, no mRNA is produced. The whole system can be expressed as:

[0121]

[0122] where G and G * represent the active and inactive states of the gene, respectively. The telegraph model explains the randomness of gene activity, including noise, changes in gene expression levels, and the dynamics of gene regulation. Figure 3 The results show that the invention can accurately solve complex chemical master equations.

[0123] Furthermore, Figure 4 shows the efficiency comparison with other methods on a more complex gene expression model according to the present invention. The present invention considers an oscillation model, which illustrates a simple genetic negative feedback loop. In this model, a gene expresses protein X. After a fixed delay τ, X undergoes a biochemical process to be converted into another protein Y. Subsequently, Y binds to the promoter of the gene and inhibits the transcription rate of X. The reaction formula of this process is as follows:

[0124]

[0125] The reaction propensities J1(Y) and J2(Y) are defined as:

[0126]

[0127] Both J1(Y) and J2(Y) use the Hill function. The synthesis rate J1(Y) depends on the concentration of the regulator S and is inhibited by protein Y. Among them, K dis the dissociation constant of the interaction between Y and the gene promoter, and p represents the Hill coefficient. For the degradation rate J2(Y), E T is the total concentration of the protease, k2 is the turnover rate, and K m represents the Michaelis constant. Figure 4 b shows the solution effect of the chemical master equation of the present invention on the oscillation model. Figure 4 c compares the sample size required for the method of the present invention and the simple multi-layer perceptron training. Figure 4 d compares the time required for the method of the present invention and the simple multi-layer perceptron training.

[0128] Figure 5 shows the extrapolation method and extrapolation effect according to the present invention for different reaction kinetic parameters; first, the present invention explores the extrapolation ability in different kinetic parameters and delay mechanisms. In the real world, this time delay τ becomes a random variable controlled by a probability distribution. In the present invention, it is assumed that τ follows a uniform distribution with an average value of In the present invention, an additional node called Attribute is introduced into the input layer of the decoder. By changing the value of Attribute, the decoder can reconstruct different probability distributions corresponding to different delay mechanisms. During training, the present invention sets Attribute = 0 to reconstruct the distribution, and sets Attribute = 1 to reconstruct the distribution. The loss function is defined as the sum of two sets of reconstruction errors and the KL divergence. A key advantage of the VAE is its ability to unravel complex data relationships within the decoder, thus enabling interpolation between data points. By utilizing the linear relationship between the uniform distribution coefficient and the attribute, the decoder can predict the probability distribution of any τ that follows a uniform distribution with an average value of Figure 5 b shows the extrapolation effect of the present invention for different reaction kinetic parameters.

[0129] Figure 6 shows the extrapolation method and extrapolation effect according to the present invention for different reaction structures; next, the present invention evaluates the extrapolation ability of different reaction topologies. Using the variational autoencoder trained in Figure 5 , the present invention directly tests its performance on the telegraph model (see Figure 6 a). Figure 6 ​b shows the prediction accuracy of the telegraph model for various kinetic parameters and delay mechanisms. A key feature of the present invention is its ability to adapt to different kinetic parameters and topologies. This adaptability indicates that the present invention can be generalized to a set of reaction systems rather than being limited to a single reaction system, such as a reaction system in which RNA molecules degrade at a fixed time. This flexibility enables the model to be applied to different biochemical scenarios, making it both robust and versatile. In addition, this ability helps to study how changes in regulatory networks or kinetic rates affect gene expression, which is crucial for understanding evolutionary processes as well as the effects of mutations or external perturbations.

[0130] Any process or method description shown in the flowchart or described otherwise herein can be understood to represent a module, segment, or portion of code including one or more executable instructions for implementing a specific logical function or process. The scope of the preferred embodiments of the present invention includes additional implementations in which functions may be executed not in the order shown or discussed, including substantially concurrently or in reverse order according to the functions involved, which should be understood by those skilled in the art to which the embodiments of the present invention pertain.

[0131] It should be understood that the various parts of the present invention can be implemented by hardware, software, firmware, or a combination thereof. In the above embodiments, multiple steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution device.

[0132] Those of ordinary skill in the art of the present technology can understand that all or part of the steps carried by the method of the above embodiments can be completed by instructing relevant hardware through a program. The program can be stored in a computer-readable storage medium. When the program is executed, it includes one or a combination of the steps of the method embodiments.

[0133] The above-mentioned storage medium can be a read-only memory, a magnetic disk, an optical disc, etc.

[0134] In the description of this specification, the descriptions referring to terms such as "one embodiment", "some embodiments", "example", "specific example", or "embodiment" etc. mean that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in a suitable manner in any one or more embodiments or examples.

[0135] Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those of ordinary skill in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention.

[0136] By adopting the method for efficiently solving the chemical master equation based on the variational autoencoder and the quantization tensor sequence of the present invention, the complexity of solving the chemical master equation of the biochemical reaction system is reduced, the sample size and training time of training are reduced, and the trained variational autoencoder can predict the dynamics of unseen time points, different reaction structures, and molecules not included in the training data, providing a powerful tool for analyzing stochastic biochemical systems.

[0137] In this specification, the present invention has been described with reference to specific embodiments thereof. However, it is obvious that various modifications and variations can still be made without departing from the spirit and scope of the present invention. Therefore, the specification and drawings should be regarded as illustrative rather than restrictive.

Claims

1. A method for efficiently solving and processing chemical master equations based on variational autoencoders and quantized tensor sequences, characterized in that, The method described above includes the following steps: S1: Using a random simulation algorithm, perform random simulation for a specific biochemical reaction to obtain the number of target molecules, and calculate the probability distribution of the target molecules as a data set; S2: Input the initial probability distribution of the target molecules into the trained variational autoencoder, and use the variational autoencoder to map the probability distribution to an equivalent propensity function, and iteratively obtain the equivalent reaction propensity function of the intermediate reaction; S3: Input different biochemical reaction kinetic parameters into the trained variational autoencoder, that is, efficiently solve the chemical master equation of the biochemical reaction, so as to obtain the probability distribution solution of the chemical master equation.

2. The method for efficiently solving the chemical master equation based on the variational autoencoder and the quantization tensor sequence according to claim 1, characterized in that, The specific steps of step S2 include: S21: Consider a reaction system containing molecule X i (i = 1, …, P) and R reactions, and classify the indices of these molecules into a set Divide the entire reaction set into two non - overlapping subsets A and B, which have R A and R B reactions respectively. Among them, the reaction kinetics in set A is the main concern, while the reaction kinetics in set B is the secondary concern; divide all molecules into three categories: where set contains the indices of molecules involved in both set A and B, contains the indices of molecules only related to set A, contains the indices of molecules only related to set B; S22: According to the partitioning result of the entire reaction system, determine the equivalent propensity function that the variational autoencoder needs to learn, and obtain a simplified chemical master equation; S23: Input the initial probability distribution of the molecules of interest into the variational autoencoder, and divide the output into two parts, which are used as the mean and variance of a Gaussian distribution; S24: Perform a reparameterization operation on the output of the variational autoencoder to obtain the input of the variational decoder, making it possible to calculate the gradient; S25: Input the reparameterized result into the decoder, so that the obtained output is added with a bias term, and finally obtain the output of the variational autoencoder, that is, the propensity function of the equivalent reaction; S26: Use the output of the variational autoencoder to calculate the probability distribution at the next time point, and repeat the above steps until the termination time of the reaction is reached; S27: After obtaining the equivalent reaction propensity function of the intermediate reaction, obtain a simplified chemical master equation, and convert the simplified chemical master equation into a quantized tensor sequence form, so as to efficiently solve the complex chemical master equation; S28: The loss function consists of the reconstruction error and the KL divergence. Input the result obtained by solving the chemical master equation back into the encoder, and use the obtained output result to calculate the KL divergence. Calculate the gap between the result obtained by solving the chemical master equation and the probability distribution obtained by the random simulation algorithm as the reconstruction error. Repeat steps S22 - S28 until the overall loss function of the network reaches the preset target, and finally obtain the trained variational autoencoder model.

3. The method for efficiently solving and processing chemical master equations based on variational autoencoders and quantization tensor sequences according to claim 1, wherein, The specific step of step S21 is: Divide the indices of all molecules, that is, P kinds of molecules, into three categories: wherein represents the complete set of all molecular indices, H is the number of molecular species involving sets A and B, M - H + 1 is the number of molecular species related only to set A, P - M + 1 is the number of molecular species related only to set B, and M, H, and P need to be set manually according to different reaction systems. The set contains the indices of the molecules involving sets A and B, contains the molecules related only to set A, contains the molecules related only to set B, then the entire reaction system is expressed as: In this reaction system, s ir and s i ir are both non-negative integers, representing the quantities of reactants and products respectively. The reaction coefficient matrix S is defined as: S ir = s ′ ir - s ir (2) To solve the curse of dimensionality and reduce the computational complexity, use a single equivalent reaction to replace the reactions in set B. This reaction is used to accurately capture their impact on the reaction kinetics in set A, and this reaction is expressed as:

4. The method for efficiently solving the chemical master equation based on the variational autoencoder and the quantization tensor sequence according to claim 3, characterized in that, The specific step of step S22 is: Assume that the equivalent tendency function of the equivalent reaction (3) is f eff , the reduced chemical master equation corresponding to the reaction system composed of the reactions in set A and the equivalent reaction (3) is given by the following formula: where n i is the number of molecules of molecule X i , S r is the r-th column of matrix S, and f r (n) is the tendency function for the r-th reaction to occur when the system is in state n, expressed as: where k r is the reaction rate of the r-th reaction, and Ω is set to 1; Utilize the universal approximation ability of the neural network and use VAE(NN θ ) to approximate f eff , where θ is the neural network parameter of VAE, and it is represented by the following formula: Assume that N represents the maximum number of all M kinds of molecules. Therefore, equation (4) is concisely expressed as: Among them, P(t) represents a vector containing all possible states, ranging from to Among them, for a single molecule, represents the probability distribution of the species, and the state transition matrix A θ (t) is defined as A θ (t) = D + N θ (t), where D contains the propensity function f of the reactions in set A r , and N θ (t) contains the effective propensity NN θ (n, t) related to the reactions in set B.

5. The method for efficiently solving and processing chemical master equations based on variational autoencoders and quantization tensor sequences according to claim 4, wherein The specific step of step S23 is: Design a suitable encoder E(·; φ), where · represents the input of the neural network and φ is the parameter of the neural network. The main task of the encoder is to map the input data to a distribution in the latent space and is used to estimate the variational distribution. Taking P(t) as the input, generate the variational distribution q(z|P(t); φ). Assume that the variational distribution q(z|P(t); φ) follows a Gaussian distribution with a diagonal covariance matrix: where I is the identity matrix, μ I and σ I are the mean and standard deviation of the Gaussian distribution, z is the latent variable, and these parameters are predicted by the encoder:

6. The method for efficiently solving the chemical master equation based on the variational autoencoder and the quantization tensor sequence according to claim 5, wherein The specific content of step S24 is as follows: Since the random variable z is sampled from the variational distribution q(z|P(t); φ), z is not deterministically related to φ. To enable the model to perform backpropagation and optimize the parameters, z is reparameterized as: z = μ I + σ I ⊙ ∈ (10) Among them, and ⊙ represent the Hadamard product. This reparameterization method converts the random relationship between z and φ into a deterministic relationship, so that the derivative of z with respect to φ can be directly calculated.

7. The method for efficiently solving and processing chemical master equations based on variational autoencoders and quantization tensor sequences according to claim 6, characterized in that The specific content of step S25 is as follows: Design a suitable decoder D(·; ψ), where · represents the input of the neural network and ψ is the parameter of the neural network. Using the reparameterized result z as the data of the decoder, obtain the output of the variational autoencoder: D(z; ψ) = p(P(t)|z; ψ) (11) To accelerate the training, a set of preset biases are added to encode prior knowledge about the reaction kinetics in set B. This bias enhances the model by providing appropriate initial parameters for the equivalent propensity function. Therefore, the final output of the VAE model is given by the following formula:

8. The method for efficiently solving and processing chemical master equations based on variational autoencoders and quantization tensor sequences according to claim 4, characterized in that, The specific content of step S27 is as follows: Since matrix A θ (t) is a high-dimensional matrix and direct representation would lead to huge storage requirements. To effectively represent this matrix, the tensor train (TT) format is used to solve the problem: Suppose there is an N-dimensional matrix A θ (t) is of size I1×I2×…×I N , and is used to represent in the format of a tensor train: A θ (i1,…,i N )=A1(i1)·A2(i2)·…·A N (i N ) (14) where each A k (i k ) is a matrix of rank r k-1 ×r k ; Suppose there is an N-dimensional matrix A θ (t) with size I1×I2×…×I N , first remap the size I of each pattern k to a power of 2, 2L k , and then reshape the tensor into a new tensor Y: Y(i1,i2,…,i N ) = A H (i1,i2,…,i N ) (15) The size of the new tensor Y is 2×2×…×2, and the index of each mode is represented in binary, so that it can be efficiently stored and solved using the QTT format.

9. The method for efficiently solving and processing chemical master equations based on variational autoencoders and quantization tensor sequences according to claim 8, characterized in that, The specific content of step S28 is as follows: The loss function of the VAE is the sum of the reconstruction error and the difference between the posterior distribution and the prior distribution in the latent space, and it can be expressed as: where β is the weight of the reconstruction error, represents the Kullback-Leibler (KL) divergence, which is used to quantify the difference between the posterior distribution q(z|P(t); φ) and the prior distribution p(z|P(t); ψ); the reconstruction error is the mean square error (MSE) between the obtained probability distribution P(t) solved and the reference distribution generated by the exact stochastic simulation algorithm: Among them, is the set of time points for collecting training data, and N shots is the total number of selected time points. The total KL divergence of all selected time points is given by the following formula: where J is the dimension of the latent space z, and are the mean and standard deviation vectors of the Gaussian distribution.

Images