Clinical decision support method based on QSP and deep reinforcement learning

Abstract

The invention discloses a clinical decision support method based on QSP and deep reinforcement learning. The method comprises the following steps: preprocessing and standardizing physiological data of a patient and existing drug parameters; constructing a QSP model of the medicine and the patient; based on the LSTM model and the ARIMA model, constructing a patient physiological state prediction model; a state space and an action space are defined, and a transition probability and a reward function are obtained; calculating a Q value and using the Q value to initialize the Q value of the DQN; accumulating experience in the experience pool by using an epsilon-greedy strategy; training the Double DQN model until the Q function is converged; and updating parameters of the Q network based on a genetic-annealing algorithm. By collecting and analyzing the detailed physiological data of the patient, a customized treatment scheme is provided for each patient, the treatment effect is improved, and unnecessary side effects are reduced.

Description

Technical Field

[0001] This application belongs to the field of clinical decision-making technology, and particularly relates to a clinical decision support method based on QSP and deep reinforcement learning. Background Art

[0002] Clinicians generally follow guidelines when prescribing medications. Due to individual differences, the treatment effects of some patients are not satisfactory. The lack of personalized treatment plans is also a dilemma in current clinical treatment. In recent years, the dose optimization plan for clinical trials has gradually shifted to using digital modeling methods to simulate the differences in individuals after receiving drug treatment in order to reduce the cost of clinical trials. For example, the QSP model quantitatively describes the pharmacokinetics (PK) and pharmacodynamics (PD) of drugs, and deep reinforcement learning selects the best plan in the environment by learning and optimizing strategies. The QSP (Quantitative Systems Pharmacology) model includes the processes of drug absorption, distribution, metabolism, and excretion (ADME) in the body, as well as the interaction between the drug and the target. The RL (Reinforcement Learning) algorithm achieves the final dose optimization by continuously optimizing decisions in a given environment.

[0003] Although the QSP model has been widely used in the process of drug research and development, relying solely on the QSP model for dose optimization still cannot achieve the best dosing plan. The LSTM (Long Short-Term Memory) algorithm has outstanding capabilities in solving non-linear problems, and the ARIMA (Autoregressive Integrated Moving Average) model has great advantages in short-term prediction. We combine the two to form a side effect index prediction model, which can comprehensively integrate the advantages of both [3]. In the field of deep reinforcement learning, the Value Iteration Network (VIN) is an important deep learning model that can be used to solve the problem of slow decision-making in Markov processes. The Double DQN model (Double Deep Q-Network) is used as the basic framework of the algorithm [4]. Compared with the traditional DQN model, it can avoid the overestimation of Q values. The genetic annealing algorithm is used to update the gradient of the neural network, and its advantages over the traditional gradient descent for updating neural network parameters lie in its global search ability and multi-objective optimization ability.

[0004] In summary, the existing drug dose optimization plans are difficult to achieve the expected effects of researchers. Inventing a deep reinforcement learning optimization algorithm based on the QSP model in the medical field will improve this dilemma, better assist in diagnosis and treatment, and help medical staff give more effective dosing plans. Summary of the Invention

[0005] This application provides a clinical decision support method based on QSP and deep reinforcement learning to solve the above technical problems.

[0006] To solve the above technical problems, a technical solution adopted in this application is: A clinical decision support method based on QSP and deep reinforcement learning, including:

[0007] Preprocess and standardize the physiological data of the patient and the existing drug parameters to obtain standard data;

[0008] Construct a QSP model of the drug and the patient; verify and optimize the QSP model as a pharmacodynamic index prediction model;

[0009] Based on the LSTM model and the ARIMA model, construct a patient physiological state prediction model, and use the patient physiological state prediction model as a side effect index prediction model;

[0010] Define the state space and action space to obtain the transition probability and reward function;

[0011] Based on the value iteration network, calculate the Q value and use the Q value to initialize the Q value of the DQN;

[0012] Use the ε-greedy strategy to accumulate experience in the experience pool, and randomly sample and output to the Double DQN model after accumulating to the set value; among them, the Double DQN model includes a Q network and a target network;

[0013] Use the data randomly sampled from the experience pool as the input of the Double DQN model, and train the Double DQN model until the Q function converges;

[0014] Update the parameters of the Q network based on the genetic-annealing algorithm.

[0015] Furthermore, the QSP pharmacodynamic index prediction model is composed of a physiological pharmacokinetic model and a pharmacodynamic model, and predicts the behavior and effects of drugs in the body by simulating the processes of drug absorption, distribution, metabolism, and excretion in the human body to generate more data that conform to the actual situation.

[0016] Furthermore, the method for constructing a patient physiological state prediction model includes:

[0017] Design the parameters in the LSTM model and randomly initialize the LSTM model;

[0018] Use the standard data as the input of the LSTM model to obtain the input gate, forget gate, output gate, and cell state of the LSTM unit; among them, each LSTM unit has a separate weight matrix and bias term to control the update of the input gate, forget gate, output gate, and cell state;

[0019] Calculate the loss function and use backpropagation to update the parameters of the LSTM model;

[0020] Use the prediction result of the LSTM model as a new feature and merge it with the original time series data into a composite vector; among them, the composite vector is used as the input of the ARIMA model.

[0021] Based on the autocorrelation function and partial autocorrelation function, determine the parameters of the ARIMA model.

[0022] Based on the parameters of the ARIMA model, construct the ARIMA model and train the ARIMA model.

[0023] Based on the trained ARIMA model, predict the standard data; combine the prediction result of the ARIMA model with the prediction result of the LSTM model.

[0024] Furthermore, based on formula (1), obtain the reward function; among them, formula (1) is:

[0025] R(s,a) = α * E(s,a) - β * A(s,a) (1);

[0026] Among them, α and β are weight coefficients used to balance the importance of drug efficacy and side effects; E(s,a) represents the drug efficacy function, which quantifies the therapeutic effect of the drug; A(s,a) represents the side effect function, which quantifies the side effects that the drug may cause.

[0027] Furthermore, a method for calculating the Q value based on the value iteration network and using the Q value to initialize the Q value of the DQN includes:

[0028] Define the discount factor γ, initialize the value function V(s), and perform value iteration update using formula (2): Among them, formula (2) is:

[0029] V(s) = max a [R(s,a) + γ * ∑ s’ P(s’|s,a) * V(s’)] (2);

[0030] Among them, V(s) is the value of state s in the next iteration step, and V(s’) is the value of state s’ in the current iteration step; max a means taking the maximum value for all possible actions a. When the value function converges, that is, when the change in the value function between two adjacent iterations is less than a certain threshold or reaches the maximum number of iterations, stop the iteration.

[0031] Based on formula (3), calculate the final Q value; among them, formula (3) is:

[0032] Q(s,a) = R(s,a) + γ * ∑ s’ P(s’|s,a) * V(s’) (3);

[0033] Among them, V(s’) is the value function of state s’ obtained through the neural network;

[0034] Set the parameters of the value iteration network, and extract the Q value corresponding to each state-action from the trained value iteration network, so that the Q value corresponding to each state-action is associated with the initial weight of the corresponding neuron in the output layer of the value iteration network, ensuring that the initialization of all weights is based on the Q value.

[0035] Furthermore, a method of using the ε-greedy strategy to accumulate experience in the experience pool and randomly sample and output to the Double DQN model after accumulating to a set value includes:

[0036] Input the initial state s, randomly generate a random number greater than or equal to 0 and less than 1, and compare the random number with the probability ε of random exploration; if the random number is less than ε, randomly select an action; if the random number is greater than or equal to ε, select the optimal action under the current estimate;

[0037] Input the state s into the Q network, obtain the Q values of each action in the state s, and select the action with the highest Q value;

[0038] Based on the QSP model and the patient physiological state prediction model, calculate the next state s’, and use the reward function to calculate the obtained reward value, and judge whether the next state s’ is the final state;

[0039] Store [s, a, r, s′, done] in the experience replay pool, where s is the current state, a is the taken action, r is the obtained reward, s′ is the next state after taking the action, done indicates whether the next state s’ is the final state, and detect the number of data in the experience replay pool; if the number of data is less than the set value, assign s′ to s, and repeat this step until the number of data is greater than the set value;

[0040] Randomly select batch_size number of [s, a, r, s′, done] from the experience pool, and the batch size determines the number of samples used for each network update.

[0041] Furthermore, a method of training the Double DQN model includes:

[0042] Based on formula (4), use the current Q evaluation network to select the action of the maximum value function; where formula (4) is:

[0043] a max = arg max a Q main (φ(s′j), a; θ) (4);

[0044] Among them, the predicted Q value output by the Q network is Q main(s, a; θ), where θ are the parameters of the Q evaluation network; φ(s) synthesizes the current state s in the current batch of data into a new vector;

[0045] Based on formulas (5)-(6), for action a max Calculate the target Q value in the target network; Formulas (5)-(6) are:

[0046] y j = r j + γQ target (φ(s′ j ), a max ; θ’) (5);

[0047] y j = r j + γQ target (φ(s′ j ), arg max a Q main (φ(s′ j ), a; θ); θ’) (6);

[0048] where r is the reward and θ’ are the parameters of the target network;

[0049] Based on formula (7), calculate the loss function:

[0050]

[0051] where L(θ) is the loss function.

[0052] Furthermore, the method for updating the parameters of the Q network based on the genetic-annealing algorithm includes:

[0053] Initialize the population, randomly generate a set of DQN parameters within a given range as chromosomes, and form an initial population with a given number of chromosomes;

[0054] Encode the parameters in the DQN into chromosomes in the form of real number coding in the genetic algorithm, and calculate the fitness of each chromosome through the fitness function; Record the chromosome with the highest fitness as the current optimal solution;

[0055] Perform the crossover operation on the chromosomes to obtain the crossover chromosomes;

[0056] Perform the mutation operation on the crossover chromosomes to obtain the mutant chromosomes;

[0057] Call the fitness function for all chromosomes including the crossover chromosomes and the mutant chromosomes, compare the chromosome with the highest fitness with the current optimal solution, and determine whether to update the current optimal solution;

[0058] Check whether the set number of iterations is reached; if so, stop the iteration and output the current optimal solution; if not, repeat the above steps;

[0059] Check whether the number of rounds of Q-network update is an integer multiple of the update frequency; if so, copy the parameters of the Q-network to the target network to achieve delayed update of the target network; if not, skip the update of the target network.

[0060] The beneficial effects of this application are as follows: By collecting and analyzing the detailed physiological data of patients, this application provides customized treatment plans for each patient, improves the treatment effect and reduces unnecessary side effects. Using the QSP pharmacodynamic index prediction model, the side effect index prediction model based on LSTM and ARIMA, and the deep reinforcement learning algorithm, the system can accurately simulate the dynamic process of drugs in patients and optimize the dosage and time of drug administration to achieve more precise treatment effects. Through automated decision support, this system helps reduce the time and effort required for doctors to formulate treatment plans, enabling doctors to focus more on the clinical treatment of patients. Reducing medical costs: Precise drug administration reduces drug waste and unnecessary medical interventions, helping to reduce the overall medical cost. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 is a schematic flowchart of an embodiment of the clinical decision support method based on QSP and deep reinforcement learning of this application;

[0062] Figure 2 is a flowchart of an embodiment of the clinical decision support method based on QSP and deep reinforcement learning of this application;

[0063] Figure 3 is Figure 1 a schematic flowchart of an embodiment of step S3 in

[0064] Figure 4 is Figure 1 a flowchart of an embodiment of step S6 in DETAILED DESCRIPTION OF THE EMBODIMENTS

[0065] To make the objectives, technical solutions and advantages of the present invention clearer, the following further describes the present invention in detail with reference to specific embodiments.

[0066] Many specific details are set forth in the following description in order to provide a thorough understanding of the present invention. However, the present invention may be implemented in other ways different from those described herein. Therefore, the present invention is not limited by the limitations of the specific embodiments disclosed in the following specification.

[0067] Refer to Figure 1-2 , Figure 1It is a schematic flowchart of an embodiment of the clinical decision-making support method based on QSP and deep reinforcement learning of the present application. The method includes:

[0068] Step S1. Preprocess and standardize the physiological data of the patient and the existing drug parameters to obtain standard data.

[0069] Specifically, collect the physiological data of the patient, the physical and chemical information of the drug, and the physiological parameters related to pharmacokinetics and pharmacodynamics of the drug, standardize and preprocess the data, so as to obtain standard data to meet the input requirements of the model.

[0070] Step S2. Construct a QSP model of the drug and the patient; verify and optimize the QSP model as a pharmacodynamic index prediction model.

[0071] Specifically, refer to Figure 3 , Figure 3 is Figure 1 a schematic flowchart of an embodiment of step S2 in

[0072] Step S201. Define the research objective and collect experimental data.

[0073] Specifically, predict the change of the blood drug concentration of aspirin in different patients over time. At the same time, collect the physical and chemical information of the drug aspirin, including lipophilicity, the fraction of the drug not bound to plasma proteins or other blood components, relative molecular mass, dissociation constant, solubility, permeability, logarithm of the oil-water partition coefficient. Collect relevant population data, including basic data such as age, gender, height, and weight. And relevant physiological parameters, including liver and kidney weights, total liver clearance rate, and renal plasma clearance rate.

[0074] Step S202. Select the model type.

[0075] Select the PBPK model in PK-sim and set the parameters of the model. Input the relevant parameters of the drug aspirin mentioned in step 1.2. Set the metabolic pathway of aspirin, add the main enzyme (CYP2C9) in the metabolic process of aspirin, and look up the initial concentration and clearance rate of CYP2C9 in the data and input them.

[0076] Step S203. Set the dosing regimen.

[0077] Specifically, determine the dose and method of input administration. The dosing interval is 24h and the number of doses is 30 times. The dosing method in this model is oral administration, single dose. Then define the dosing interval and number of times. Define the characteristics of individuals in the selected age range in the individual section, including age, weight, height, and gender ratio.

[0078] Step S204. Create a new simulation, add the section set in the above steps to the simulation, and PK-sim will simulate the pharmacokinetic process of the drug in the body according to the input parameters and dosing regimen.

[0079] Step S205. Analyze the results. After the simulation ends, the concentration-time curve of the drug in the body can be obtained.

[0080] Step S206. Validate the model. Compare the data obtained from experimental observations or found in the literature with the data predicted by the model, and verify the model by comparing the consistency between the observed data (Observed Data) and the predicted data (Predicted Data).

[0081] Step S207. Optimize the model. Consult the literature to find the parameters and physiological indicators that affect or are related to the pharmacokinetic process of aspirin in the body, including lipophilicity, permeability, total hepatic clearance, renal plasma clearance, dose, and body weight. Add them to the "Parameter Identification" section for parameter optimization, and then "Transfer to Simulation" the optimized values and apply them to the model.

[0082] Step S208. Construct and validate the population model. Create a population based on the age range of the corresponding individuals in the population section, construct a population PBPK model, and verify the absorption and distribution of the drug in this population.

[0083] Step S209. Analyze the results. Analyze the C_max, AUC, t_max, and half-life of the optimized model, and compare the consistency between the observed data and the predicted data.

[0084] Step S210. Sensitivity analysis. Create a sensitivity analysis, add all the parameters to the sensitivity analysis to determine which parameters have a greater impact on the model output.

[0085] Step S211. Iterative optimization. Continuously calibrate and optimize the model by inputting a large amount of clinical data or experimental data to make the model more accurate.

[0086] Step S212. Define the research objective. Predict the changes in platelet aggregation rate in the body after different patients take a small dose of aspirin orally. According to the mutual relationship between platelet aggregation rate, TXA2, and COX-1 concentration, quantitatively convert the required change in platelet aggregation rate into a change in COX-1 concentration, and then obtain the change in COX-1 concentration through PK and PD models. The quantitative conversion is based on the formula:

[0087] Agg = β * α * COX-1;

[0088] where α and β are proportionality constants obtained from experimental data.

[0089] Step S213. Export the PK model built in PK-sim in.json format, import it into MATLAB, and convert the read binary data into a character vector to obtain a data structure that MATLAB can process. Extract the relevant parameter values required from the processed data. The parsing of the JSON string is read and parsed using the jsondecode function in MATLAB.

[0090] Step S214. Select a suitable PD model according to the way and result of the pharmacodynamic effect of aspirin in the body. Here, the indirectly connected S-shaped E max model is selected. The relevant parameter data includes the administered dose of aspirin, blood drug concentration, and time data.

[0091] Step S215. Define the E max model by adding the extracted blood drug concentration values to the corresponding effect values, and set initial guess values for the relevant parameters in the model. The relevant parameters include E min (effect baseline without drug), E max (theoretical maximum of drug effect), EC 50 (concentration when the drug concentration reaches half of the maximum effect) and Hill coefficient.

[0092] Step S216. Set the boundaries of the parameters to ensure that the parameter values are within a reasonable range during the fitting process.

[0093] Step S217. Use the lsqcurvefit function in MATLAB for model fitting, minimize the squared difference between the model predicted values and the actual observed values, and output the fitted parameter values.

[0094] Step S218. Use the parameter values output in Step S217, plus the blood drug concentration values extracted in Step S213, to establish an S-shaped Emax model according to the equation. The equation used is:

[0095]

[0096] where n is the Hill coefficient. Obtain the change of COX-1 activity with blood drug concentration.

[0097] Step S219. Combine the law of the change of COX-1 activity with time, and convert it to obtain the curve of the change of COX-1 activity with time. Plot the curve of the change of COX-1 activity with time, set the X-axis and Y-axis labels and the image title, enable the grid for easy observation. Calculate the peak and minimum values within each dosing period, mark them in the image and display the peak and minimum values.

[0098] Step S3. Based on the LSTM model and the ARIMA model, construct a patient physiological state prediction model, and use the patient physiological state prediction model as the side effect index prediction model.

[0099] Specifically, step S3 includes:

[0100] S31. Design the parameters in the LSTM model and randomly initialize the LSTM model.

[0101] Specifically, collect the patient's basic physiological data such as gender, age, height, weight, and parameters such as blood glucose and blood lipids, input them into the LSTM model, and randomly initialize the LSTM model.

[0102] S32. Use the standard data as the input of the LSTM model to obtain the input gate, forget gate, output gate, and cell state of the LSTM unit; among them, each LSTM unit has its own weight matrix and bias term, which are used to control the update of the input gate, forget gate, output gate, and cell state.

[0103] Specifically, the number of layers of the LSTM is designed to be 2 layers, the number of units in the first layer of the LSTM is 128, and the number of units in the second layer of the LSTM is 64. For the return sequence (whether the LSTM layer outputs its internal state at each time step), the first layer is designed to be True, and the second layer is designed to be False. Then, use the init_lstm_state function to initialize the hidden state and cell state of the LSTM. The learning rate is designed to be 0.001, which can be adjusted during the training process. The Adam optimizer is used, which is an adaptive learning rate optimization algorithm that can adjust the learning rate according to the dynamic changes during the training process. The discount factor is set to γ = 0.99. The discount factor γ determines the discount degree of future rewards and affects the agent's emphasis on immediate and future rewards. The number N of batches processed at one time is designed to be 10.

[0104] The core of the LSTM unit lies in its three gating mechanisms (input gate, forget gate, and output gate) and the cell state. Each LSTM unit has its own weight matrix and bias term, which are used to control the update of these gates and the cell state.

[0105] The data processing flow is as follows:

[0106] Input gate:

[0107] The input data at the current time step: x t The hidden state at the previous time step: h t-1 The purpose of the input gate is to determine how much new information enters the cell state. The output of the input gate is a vector i obtained through a sigmoid activation function. t The calculation formula of the input gate is:

[0108] i t = σ(W i [x t ,h t-1 +b i )

[0109] where W i is the weight matrix, b i is the bias term, and σ i is the sigmoid activation function.

[0110] Forget gate:

[0111] The purpose of the forget gate is to determine how much of the cell state from the previous time step to retain. The output of the forget gate is a vector f t obtained through a sigmoid activation function. The calculation formula for the forget gate is:

[0112] f t = σ(W f [x t ,h t-1 +b f )

[0113] where W f is the weight matrix and b f is the bias term.

[0114] Calculate the new candidate cell state, obtaining a vector in the range [-1, 1] using the tanh activation function. The calculation formula for the new information is:

[0115]

[0116] where W c is the weight matrix, b c is the bias term, and tanh is the hyperbolic tangent activation function.

[0117] Cell state update:

[0118] When updating the cell state c t the output f t of the forget gate and the output i t of the input gate as well as the new information are considered. The cell state update formula is:

[0119]

[0120] where ⊙ represents element-wise multiplication.

[0121] Output gate:

[0122] The purpose of the output gate is to determine how much information to output from the cell state. The output of the output gate is a vector o obtained through a sigmoid activation function. t . The calculation formula of the output gate is:

[0123] o t = σ(W o (x t ,h t-1 +b o );

[0124] Among them, W o is the weight matrix, and b o is the bias term.

[0125] Hidden state calculation:

[0126] The hidden state h t is determined by the output o of the output gate t and the cell state c passed through the tanh activation function t . The calculation formula of the hidden state is:

[0127] h t = o t ■tanh(c t );

[0128] Among them, this is the final output of the LSTM unit.

[0129] S33. Calculate the loss function and use backpropagation to update the parameters of the LSTM model.

[0130] Specifically, use the output of the network and the actual target value to calculate the loss function. The formula of the loss function is:

[0131]

[0132] Among them, N represents the number of samples in the batch, X and Y respectively represent the actually measured blood glucose and blood lipid levels, respectively represent the blood glucose and blood lipid levels predicted by the model. And use this loss function to calculate the gradients of the network output, the output gate, and the hidden state in sequence, and backpropagate the gradient of the hidden state to each time step to calculate the gradients of the input gate, the forget gate, the cell state, and the input at each time step, and finally realize parameter update.

[0133] S34. Use the prediction result of the LSTM model as a new feature and merge it with the original time series data into a composite vector; among them, the composite vector is used as the input of the ARIMA model.

[0134] Specifically, the prediction result of the LSTM model is used as a new feature and merged with the original time series data to form a composite vector α(s) as the input of the ARIMA model.

[0135] S35. Based on the autocorrelation function (ACF) and partial autocorrelation function (PACF), determine the parameters of the ARIMA model.

[0136] ACF (Autocorrelation Function): Calculate the mean of the time series data for subsequent calculations. For each time point, calculate its lag values with itself and other time points. Calculate the Pearson correlation coefficient between the value at the current time point and all lag values. Plot the calculated correlation coefficients on a graph, with the horizontal axis representing the lag order and the vertical axis representing the correlation coefficient. Determine the significance of the ACF to be set at 0.05. If the correlation coefficient is greater than the significance level at a certain lag order, it is considered that there is a significant correlation between the sequence and its own lag values. Through ACF analysis, the parameters of the AR (Autoregressive) part in the ARIMA model can be determined. ACF shows the correlation between the sequence and its own lag values, which helps to identify which lag values have a significant impact on the current value.

[0137] PACF (Partial Autocorrelation Function): First, calculate the ACF of the time series. For each time point, calculate its lag values with itself and other time points. Calculate the Pearson correlation coefficient between the value at the current time point and all lag values. Use the ACF values to eliminate the correlation between the sequence and its lag values. Add the product of the ACF values of the current time point and all previous lag values. Plot the calculated correlation coefficients on a graph, with the horizontal axis representing the lag order and the vertical axis representing the correlation coefficient. Similar to the ACF graph, the significance level of the PACF graph is also set at 0.05. If the correlation coefficient is greater than the significance level at a certain lag order, it is considered that there is a significant correlation between the sequence and its own lag values. PACF helps to determine the parameters of the MA (Moving Average) part. By eliminating the correlation between the sequence and its lag values, PACF helps to identify which lag values' influence on the current value has disappeared.

[0138] Determine the parameters (p, d, q) of the ARIMA model: Based on the analysis of ACF and PACF, we can determine the parameters of the ARIMA model. Among them, p is the order of the AR part, d is the order of differencing, and q is the order of the MA part.

[0139] S36. Based on the parameters of the ARIMA model, construct the ARIMA model and train the ARIMA model;

[0140] Specifically, based on the identified parameters (p, d, q) of the ARIMA model, an ARIMA model is constructed. The statsmodels library in Python is used to estimate the model parameters. Check whether the model residuals are white noise. If the residuals are white noise, it indicates that the model can fit the data well without remaining periodicity or trend. In this application, the Ljung-Box test in statistical tests is used to verify whether the residuals are white noise.

[0141] S37. Based on the trained ARIMA model, predict the standard data; combine the prediction results of the ARIMA model with the prediction results of the LSTM model.

[0142] Specifically, use the trained ARIMA model to predict the integrated data. Combine the prediction results of the ARIMA model with the prediction results of the LSTM model. This can combine the advantages of the two models. The ARIMA model processes linear trends and seasonality, while the LSTM model processes non-linearity and long-term dependencies.

[0143] Step S4. Based on the defined state space and action space, obtain the transition probability and the reward function.

[0144] Specifically, define the state space (S), the action space (A), design the transition probability (P) and the reward function; where the reward function is the following formula (1):

[0145] R(s,a) = α * E(s,a) - β * A(s,a) (1);

[0146] Among them, α and β are weight coefficients used to balance the importance of drug efficacy and side effects; E(s,a) represents the drug efficacy function, quantifying the therapeutic effect of the drug; A(s,a) represents the side effect function, quantifying the side effects that the drug may cause.

[0147] Step S5. Based on the value iteration network, calculate the Q value and use the Q value to initialize the Q value of the DQN.

[0148] Specifically, step S5 includes:

[0149] Step S51. Collect the basic data of the patient, including the patient's height, weight, age, gender, blood sugar, platelet aggregation rate, and blood lipid data. Normalize the values of the patient's height, weight, age, gender, blood sugar, platelet aggregation rate, and blood lipid indicators to the range [0,1] using the maximum-minimum normalization method, perform clustering using the hierarchical clustering algorithm, determine the optimal number of clusters through the silhouette coefficient and put it into the state space; record the different doses of each drug and put them into the action space.

[0150] Step S52. Design the transition probability (P(s′∣s,a)) based on the state space and the action space. The reward function formula is as follows:

[0151]

[0152] where PAR represents the platelet aggregation rate, NPAR represents the baseline platelet aggregation rate, B t and B normal represent the blood lipid concentrations in the current state and the normal state respectively, G t and G normal represent the blood glucose concentrations in the current state and the normal state respectively. α is a parameter that adjusts the impact of the change in platelet aggregation rate on the reward, and β and λ are parameters that adjust the impacts of the side - effect scores of blood lipid and blood glucose on the reward.

[0153] Step S53. Randomly initialize the value iteration network (VIN), set the discount factor γ = 0.9, and perform value iteration update using formula (2). Among them, formula (2) is:

[0154] V(s) = max a [R(s,a) + γ * ∑ s, P(s’|s,a) * V(s’)] (2);

[0155] where V(s) is the value of state s in the next iteration step, and V(s′) is the value of state s′ in the current iteration step; max a represents taking the maximum value for all possible actions a. When the value function converges, that is, when the change in the value function between two adjacent iterations is less than a certain threshold or reaches the maximum number of iterations, stop the iteration.

[0156] Step S54. Calculate the final Q value based on formula (3). Among them, formula (3) is:

[0157] Q(s,a) = R(s,a) + γ * ∑ s, P(s’|s,a) * V(s’) (3);

[0158] where V(s′) is the value function of state s′ obtained through the neural network.

[0159] Step S55. Set the parameters of the value iteration network, extract the Q value corresponding to each state - action from the trained value iteration network, so that the Q value corresponding to each state - action is associated with the initial weight of the corresponding neuron in the output layer of the value iteration network, and ensure that the initialization of all weights is based on the Q value.

[0160] Specifically, set the parameters required for the operation of the value iteration network, including: learning rate α = 0.1, discount factor γ = 0.99, exploration rate ε = 40%, and target network parameter update frequency C = 5. Redefine the state space as the patient's platelet aggregation rate, blood glucose concentration, and blood lipid concentration.

[0161] Extract the Q-values of each state-action pair from the trained value iteration network (VIN). These Q-values are part of the optimal policy learned by the value iteration network (VIN). Associate the Q-values of each state-action pair with the initial weights of the corresponding neurons in the output layer of the network, ensuring that the initialization of all weights is based on the Q-values rather than random initialization.

[0162] Step S6. Use the ε-greedy strategy to accumulate experiences in the experience pool. After accumulating to the set value, randomly sample and output to the Double DQN model; among them, the Double DQN model includes a Q-network and a target network.

[0163] Specifically, refer to Figure 4 , Step S6 includes:

[0164] Step S61. Input the initial state s, randomly generate a random number greater than or equal to 0 and less than 1, and compare the random number with the random exploration probability ε; if the random number is less than ε, randomly select an action; if the random number is greater than or equal to ε, select the optimal action under the current estimate.

[0165] Specifically, input an initial state s, randomly generate a random number between 0 (including 0) and 1 (excluding 1), and compare it with ε (representing the probability of random exploration). If the random number is less than ε, the agent explores, that is, randomly selects an action. If the random number is greater than or equal to ε, the agent exploits, that is, selects the optimal action under the current estimate. Input the state s into the Q-network to obtain the Q-values of each action in this state, and then select the action with the highest Q-value. Call the QSP pharmacodynamic index prediction model and the side effect index prediction model based on LSTM and ARIMA to calculate the next state s', and calculate the obtained reward using the reward function.

[0166] Step S62. Input the state s into the Q-network to obtain the Q-values of each action in the state s, and select the action with the highest Q-value.

[0167] Step S63. Based on the QSP model and the patient physiological state prediction model, calculate the next state s', and calculate the obtained reward value using the reward function. Determine whether the next state s' is the final state;

[0168] Step S64. Store [s, a, r, s′, done] into the experience replay pool, where s is the current state, a is the action taken, r is the obtained reward, s′ is the next state after taking the action, and done indicates whether the next state s′ is the final state. Detect the number of data in the experience replay pool; if the number of data is less than the set value, assign s′ to s and repeat this step until the number of data is greater than the set value.

[0169] Step S65. Randomly select batch_size number of [s, a, r, s′, done] from the experience pool. The batch size determines the number of samples used for each network update.

[0170] Among them, the batch size is set to batch_size = 8. The batch size determines the number of samples used for each network update. And synthesize the current state s in this batch of data into a new vector φ(s).

[0171] The currently constructed Q-network is the evaluation Q-network, which is used to select actions and estimate Q-values. Completely copy the structure and parameters of the Q-network to generate a target network, which is used to generate target Q-values to stabilize the training process.

[0172] Step S7. Use the randomly sampled data in the experience pool as the input of the Double DQN model and train the Double DQN model until the Q-function converges;

[0173] Specifically, for the inherent defect of Q-Learning itself - overestimation, we use the Double DQN algorithm to calculate the Q-value to eliminate the bias. The predicted Q-value output by the Q-network is denoted as Q main (s, a; θ), where θ represents the parameters of the Q evaluation network.

[0174] Based on formula (4), use the current Q evaluation network to select the action of the maximum value function; where formula (4) is:

[0175] a max = arg max a Q main (φ(s′ j ), a; θ) (4);

[0176] Among them, the predicted Q-value output by the Q-network is Q main (s, a; θ), θ are the parameters of the Q evaluation network; φ(s) is the current state s in the current batch of data synthesized into a new vector;

[0177] Based on formulas (5)-(6), for the action a maxCalculate the target Q value in the target network; the formulas (5)-(6) are as follows:

[0178] y j = r j + γQ target (φ(s′ j ), a max ; θ’) (5);

[0179] y j = r j + γQ target (φ(s′ j ), arg max a Q main (φ(s′ j ), a; θ); θ’) (6);

[0180] where r is the reward and θ’ is the parameter of the target network;

[0181] Based on formula (7), calculate the loss function:

[0182]

[0183] where L(θ) is the loss function.

[0184] Update the parameters θ of the current network by updating all the parameters of the Q network through the genetic simulated annealing algorithm, so that the predicted Q value is close to the target Q value.

[0185] Step S8. Update the parameters of the Q network based on the genetic-annealing algorithm.

[0186] Specifically, step S8 includes:

[0187] Step S81. Initialize the population, randomly generate a set of DQN parameters within a given range as chromosomes, and form an initial population with a given number of chromosomes.

[0188] Specifically, initialize the population, randomly generate a set of DQN parameters within a given range as chromosomes, and form an initial population with a given number of chromosomes. Set the range of weights to (-1, 1), the range of biases to (-1, 1), and the number of individuals in the population to 200. Set an initial temperature of 1000 for the simulated annealing algorithm, and the attenuation rate K of the object annealing temperature is 0.99. The annealing strategy of the object temperature is T(n + 1) = K·T(n), where n is the number of iterations.

[0189] For each individual in the population, evaluate its fitness. The fitness function formula is:

[0190] Fitness = (θ) = -L(θ);

[0191] Among them, \(L(\theta)\) represents the mean square error loss function.

[0192] Metropolis criterion formula: The probability of accepting a new solution is determined by the following formula:

[0193]

[0194] Among them, \(E(n)\) represents the energy of the solution in state \(n\) (i.e., the value of the fitness function), and \(T\) is the current temperature.

[0195] Step S82. Encode the parameters in DQN into chromosomes in the genetic algorithm in the form of real number coding, and calculate the fitness of each chromosome through the fitness function; record the chromosome with the highest fitness as the current optimal solution.

[0196] Specifically, encode the parameters in DQN into chromosomes in the genetic algorithm in the form of real number coding, call the fitness function to calculate the fitness of each chromosome, find the chromosome with the highest fitness and record it as the current optimal solution, and use the roulette wheel selection method to generate the first generation. The formula for the selection probability is:

[0197]

[0198] where \(f\) i represents the fitness value of individual \(i\).

[0199] Step S83. Perform a crossover operation on the chromosomes to obtain crossover chromosomes.

[0200] Specifically, perform the crossover operation. Use the uniform crossover method for this operation. The genes on each genome of two paired individuals are exchanged with the same crossover probability, thus forming two new individuals. The crossover probability is set to 40%. After the crossover is completed, call the Metropolis criterion formula for the chromosomes that have undergone crossover to determine whether the solution can be accepted.

[0201] Step S84. Perform a mutation operation on the crossover chromosomes to obtain mutant chromosomes.

[0202] Specifically, perform the mutation operation. Use uniform mutation. Generate a random number uniformly distributed within a certain interval for each gene value of the individual, and replace the original gene value with this random number with a certain probability. The mutation probability is set to 1%. After the mutation is completed, call the Metropolis criterion formula for the chromosomes that have undergone mutation to determine whether the solution can be accepted.

[0203] Step S85. Call the fitness function for all chromosomes including the crossover chromosomes and the mutated chromosomes, obtain the chromosome with the highest fitness, compare it with the current optimal solution, and determine whether to update the current optimal solution.

[0204] Step S86. Detect whether the set number of iterations is reached; if so, stop the iteration and output the current optimal solution; if not, repeat the above steps.

[0205] Step S87. Detect whether the number of rounds of Q-network update is an integer multiple of the update frequency C; if so, copy the parameters of the Q-network to the target network to achieve delayed update of the target network; if not, skip the update of the target network.

[0206] The beneficial effects of the present invention are as follows:

[0207] 1. By collecting and analyzing the detailed physiological data of patients, the system can provide a customized treatment plan for each patient, thereby improving the treatment effect and reducing unnecessary side effects. Using the QSP pharmacodynamic index prediction model, the side effect index prediction model based on LSTM and ARIMA, and the deep reinforcement learning algorithm, the system can accurately simulate the dynamic process of drugs in patients and optimize the dosage and time of drug administration to achieve a more accurate treatment effect.

[0208] 2. Through automated decision support, the system helps to reduce the time and effort required for doctors to formulate treatment plans, enabling doctors to focus more on the clinical treatment of patients. Reducing medical costs: Precise drug administration reduces drug waste and unnecessary medical interventions, helping to reduce the overall medical cost.

[0209] 3. The present application uses a side effect index prediction model based on LSTM and ARIMA, which retains both the ability of the LSTM algorithm to solve non-linear problems and combines the advantages of the ARIMA model in short-term prediction.

[0210] 4. By using the value iteration network (VIN) to calculate Q-values and initialize Double DQN, the present application reduces the consumption of computing resources and accelerates the training of Double DQN. Using the Double DQN model as the basic framework of the algorithm, its experience replay technology helps to improve the utilization efficiency of data. By learning from randomly sampled past experiences, the model can better generalize to new data. Compared with the traditional DQN model, it avoids the occurrence of overestimation of Q-values.

[0211] 5. The present application updates the gradient of the neural network through a genetic annealing algorithm, which has outstanding advantages over the traditional gradient descent for updating the neural network parameters in terms of global search ability, adaptability, derivative-free requirement, balance between exploration and exploitation, multi-objective optimization ability, and parallel computing ability.

[0212] The above are only the embodiments of the present application, and do not thus limit the patent scope of the present application. Any equivalent structure or equivalent process transformation made by using the content of the specification and drawings of the present application, or directly or indirectly applied in other related technical fields, shall similarly be included within the patent protection scope of the present application.

Claims

1. A clinical decision support method based on QSP and deep reinforcement learning, characterized in that: include Preprocess and standardize the patient's physiological data and existing drug parameters to obtain standard data; Constructing a QSP model of drugs and patients; validating and optimizing the QSP model as a drug efficacy index prediction model; Based on the LSTM model and the ARIMA model, a patient physiological state prediction model is constructed, and the patient physiological state prediction model is used as a side effect indicator prediction model; Define the state space and action space, obtain the transition probability and reward function; Based on the value iteration network, the Q value is calculated and used to initialize the Q value of the DQN; Using the ε-greedy strategy to accumulate experience in the experience pool, after accumulating to a set value, randomly sampling and outputting it to the Double DQN model; wherein the Double DQN model includes a Q network and a target network; Using the randomly sampled data in the experience pool as input to the Double DQN model, and training the Double DQN model until the Q function converges; Based on the genetic-annealing algorithm, the parameters of the Q network are updated.

2. The method according to claim 1, characterized in that The QSP drug efficacy index prediction model is composed of a physiological pharmacokinetic model and a pharmacodynamic model. It predicts the behavior and effect of drugs in the body by simulating the absorption, distribution, metabolism and excretion of drugs in the human body to generate more data that conforms to the actual situation.

3. The method according to claim 1, characterized in that A method for constructing a patient physiological state prediction model, comprising: Design parameters in an LSTM model and randomly initialize the LSTM model; The standard data is used as the input of the LSTM model to obtain the input gate, forget gate, output gate and cell state of the LSTM unit; wherein each LSTM unit has a separate weight matrix and bias term for controlling the update of the input gate, the forget gate, the output gate and the cell state; Calculate the loss function and update the parameters of the LSTM model using backpropagation; The prediction result of the LSTM model is used as a new feature and combined with the original time series data into a composite vector; wherein the composite vector is used as the input of the ARIMA model; Determining the parameters of the ARIMA model based on the autocorrelation function and the partial autocorrelation function; Based on the parameters of the ARIMA model, construct the ARIMA model and train the ARIMA model; Based on the trained ARIMA model, the standard data is predicted; and the prediction result of the ARIMA model is combined with the prediction result of the LSTM model.

4. The method according to claim 1, characterized in that: Based on formula (1), the reward function is obtained; wherein the formula (1) is: R(s,a)=α*E(s,a)-β*A(s,a) (1); Among them, α and β are weight coefficients used to balance the importance of efficacy and side effects; E(s,a) represents the efficacy function, which quantifies the therapeutic effect of the drug; A(s,a) represents the side effect function, which quantifies the side effects that the drug may cause.

5. The method according to claim 1, characterized in that The method of calculating a Q value based on a value iteration network and using the Q value to initialize a Q value of a DQN comprises: Define the discount factor γ, initialize the value function V(s), and perform value iteration update using formula (2): Wherein, formula (2) is: V(s)=max a [R(s,a)+γ*∑ s' P(s'|s,a)*V(s')] (2); Where V(s) is the value of state s in the next iteration step, V(s') is the value of state s' in the current iteration step; max a It means taking the maximum value of all possible actions a. When the value function converges, that is, when the change of the value function between two adjacent iterations is less than a certain threshold or reaches the maximum number of iterations, the iteration stops. Based on formula (3), the final Q value is calculated; wherein the formula (3) is: Q(s,a)=R(s,a)+γ*∑ s' P(s'|s,a)*V(s') (3); Among them, V(s′) is the value function of state s′ obtained by the neural network; The parameters of the value iteration network are set, and the Q value corresponding to each state-action is extracted from the trained value iteration network, so that the Q value corresponding to each state-action is associated with the initial weight of the neuron corresponding to the output layer of the value iteration network, ensuring that the initialization of all weights is based on the Q value.

6. The method according to claim 1, characterized in that The method of using the ε-greedy strategy to accumulate experience in the experience pool, and randomly sampling and outputting it to the Double DQN model after accumulating to a set value, includes: Input the initial state s, randomly generate a random number greater than or equal to 0 and less than 1, and compare the random number with the probability of random exploration ε; if the random number is less than ε, randomly select an action; if the random number is greater than or equal to ε, select the optimal action under the current estimate; Input the state s into the Q network, obtain the Q value of each action under the state s, and select the action with the highest Q value; Based on the QSP model and the patient physiological state prediction model, the next state s' is calculated, and the reward value calculated by the reward function is used to determine whether the next state s' is the final state; Store [s, a, r, s′, done] into the experience replay pool, where s is the current state, a is the action taken, r is the reward obtained, s′ is the next state after taking the action, and done indicates whether the next state s′ is the final state. Check the number of data in the experience replay pool; if the number of data is less than the set value, assign s′ to s, and repeat this step until the number of data is greater than the set value; Randomly select batch_size number of [s,a,r,s′,done] from the experience pool. The batch size determines the number of samples used each time the network is updated.

7. The method according to claim 1, characterized in that The method of training the Double DQN model comprises: Based on formula (4), the action of the maximum function is selected using the current Q evaluation network; wherein the formula (4) is: a max =arg max a Q main (φ(s′ j ),a;i) (4); Among them, the predicted Q value output by the Q network is Q main (s, a; θ), θ is the parameter of the Q evaluation network; φ(s) is synthesized into a new vector for the current state s in the current batch of data; Based on formulas (5)-(6), in action a max The target Q value is calculated in the target network; the formulas (5)-(6) are: y j =r j +γQ target (φ(s′ j ),a max ;θ′) (5); y j =r j +γQ target (φ(s′ j ),arg max a Q main (φ(s′ j ),a;θ);θ') (6); Where r is the reward, θ' is the parameter of the target network; Based on formula (7), the loss function is calculated: Among them, L(θ) is the loss function.

8. The method according to claim 1, characterized in that The method for updating the parameters of the Q network based on the genetic-annealing algorithm comprises: Initialize the population, randomly generate a set of DQN parameters as chromosomes within a given range, and form the initial population with a given number of chromosomes; The parameters in the DQN are encoded as chromosomes in the genetic algorithm in the form of real number encoding, and the fitness of each chromosome is calculated by a fitness function; the chromosome with the highest fitness is recorded as the current optimal solution; Performing a crossover operation on the chromosome to obtain a crossover chromosome; Performing a mutation operation on the crossover chromosome to obtain a mutated chromosome; Calling the fitness function for all chromosomes including the crossover chromosome and the mutated chromosome, obtaining the chromosome with the highest fitness and comparing it with the current optimal solution, and determining whether the current optimal solution needs to be updated; Check whether the set number of iterations has been reached; if so, stop the iteration and output the current optimal solution; if not, repeat the above steps; Detect whether the number of rounds of the Q network update is an integer multiple of the update frequency; if so, copy the parameters of the Q network to the target network to achieve delayed update of the target network; if not, skip the update of the target network.

Images