A markov chain monte carlo sampling method based on double constraint set projection mechanism

Abstract

This invention discloses a Markov chain Monte Carlo sampling method based on a dual-constraint set projection mechanism, which solves the problems of poor target distribution adaptability, slow convergence, and sensitivity to initial values ​​in traditional MCMC sampling in personalized drug delivery scenarios. The method includes: sampling initialization, using the ppk+MAPB algorithm to construct a four-layer nested target probability distribution model; determining a multivariate normal distribution as the proposed distribution; introducing a dual-constraint mechanism of Euclidean distance constraint set and orientation constraint set during iterative sampling, mapping candidate samples to the intersection of the two sets through a projection operator, and sampling numerical parameters θ and covariance matrix Σ in stages using the MH criterion; and performing convergence judgment. This invention improves the clinical effectiveness and accuracy of the sampled samples, achieves faster convergence, and can be widely applied to parameter estimation and protocol optimization scenarios in personalized drug delivery.

Description

Technical Field

[0001] This invention relates to the fields of probability statistics, machine learning, and personalized drug delivery, specifically to a Markov chain Monte Carlo sampling method based on a double-constraint set projection mechanism, applicable to a variety of complex fields. Background Technology

[0002] Personalized drug delivery is a core branch of precision medicine. Its core principle is to optimize drug dosage and administration regimens based on individual patient characteristics, achieving a "one person, one plan" approach to reduce adverse drug reactions and improve treatment outcomes. Currently, MCMC algorithms are frequently applied to personalized drug delivery; however, traditional MCMC sampling has significant limitations in its application to personalized drug delivery scenarios.

[0003] First, the construction of target distributions lacks specificity, often employing a single distribution model, which cannot adapt to the high-dimensional and complex situations in the actual scenario of individualized drug administration.

[0004] Second, traditional MCMC defines the target distribution only through likelihood × prior, without clear clinical parameter boundary constraints. It is easy to sample abnormal parameters that do not conform to physiological logic. Although these parameters meet the probability distribution, they have no clinical significance and will lead to errors in subsequent dosing regimen calculations.

[0005] Third, in traditional MCMC sampling, the iteration direction is unconstrained when generating candidate samples through random walks, which easily leads to repeated oscillations, resulting in an increase in the number of samplings, a slower convergence speed, and a higher proportion of invalid samples, thus increasing the time consumption.

[0006] Fourth, traditional MCMC is sensitive to initial values. If the initial values ​​deviate from the target distribution, it is easy to get trapped in local optima, resulting in large differences in the posterior distribution obtained with different initial values, which affects the stability of the dosing regimen.

[0007] Fifth, there are random errors in clinical TDM testing. Traditional MCMC may sample incorrect parameters that fit outliers because the likelihood function is sensitive to outliers.

[0008] Therefore, while retaining the advantages of the MCMC framework, this invention introduces a dual-constraint projection mechanism to address the specific needs of personalized drug delivery scenarios, thereby overcoming the limitations of traditional MCMC sampling in this field, optimizing the target distribution construction and sample screening mechanism, and improving the accuracy and reliability of personalized drug delivery parameter sampling. Summary of the Invention

[0009] The purpose of this invention is to provide an improved sampling method to address the core needs of personalized drug delivery scenarios. Based on known population and pharmacokinetic model information, and combined with the multi-level parameter relationships in personalized drug delivery, the method optimizes the target distribution construction and sample screening mechanism, achieves directional shrinkage of the sampling space, and provides precise design for personalized drug delivery protocols.

[0010] Technical solution

[0011] To achieve the above objectives, the present invention employs a sampling method suitable for personalized drug administration, characterized by comprising the following steps:

[0012] S1: Sampling initialization, determine the core sampling parameters in the individualized dosing scenario, construct the target probability distribution model using the hierarchical analysis method of population pharmacokinetic and maximum a posteriori Bayesian estimation (ppk+MAPB) algorithm, and initialize the initial state of the Markov chain;

[0013] S2: Determine the proposed distribution and select an appropriate proposed distribution and proposed distribution parameters based on the distribution characteristics of individualized dosing parameters;

[0014] S3: Markov chain iterative sampling, based on the Metropolis-Hastings criterion, combined with the proposal distribution, generates sampling samples, adds dual set filtering and projection operations to optimize samples, and records the chain state sequence;

[0015] S4: Convergence judgment. Use conventional convergence judgment indicators to monitor the stability of the chain, determine the number of convergence iterations, and stop invalid iterations.

[0016] S5: Sample screening, performing autocorrelation removal on the converged samples to obtain effective samples that conform to the target distribution;

[0017] S6: Sampling result output and verification. Output the valid samples and statistical characteristics, verify the consistency between the samples and the target distribution, provide data support for the optimization of individualized dosing regimens, and complete the entire sampling process.

[0018] Furthermore, step S1 specifically includes:

[0019] S11: Obtain the user's personalized drug delivery sampling requirements, including the target distribution type, parameter dimension d, target distribution probability density function f(x), total number of samples N, and initial iteration threshold. Maximum iteration threshold The target distribution is constructed using the hierarchical analysis method of the ppk+MAPB algorithm, building a multi-layered nested target probability distribution model. The specific hierarchical structure is as follows:

[0020] 1. First layer (observation layer): Let... obj is the observed value of individual i at time j, specifically the observed blood drug concentration in a personalized drug administration scenario. pre is the model-predicted blood drug concentration at time j for individual i. It is a random error term and It follows a normal distribution with a mean of 0 and a variance of σ², based on individual random effects. = The relationship between observed and predicted values ​​is obtained:

[0021] ln obj=ln pre+ ,

[0022] That is, ln obj follows a mean of ln pre, a normal distribution with variance σ²; where ln pre is determined by the individual's core pharmacokinetic parameters, namely clearance rate CL and volume of distribution V, i.e., ln pre=h( h(·) is a preset concentration calculation function, defined as:

[0023] h(·)=

[0024] in It can be represented as follows:

[0025]

[0026] 2. Second layer (individual parameter layer): Based on the population pharmacokinetic parameter assumptions, the i-th patient is taken as a sample in the population, and the core pharmacokinetic parameter vector is ( The transpose of ) follows a two-dimensional normal distribution, i.e.

[0027]

[0028] ( The transpose of ) follows a two-dimensional normal distribution with mean μ and covariance matrix Σ, where μ is the overall mean of the population pharmacokinetic parameters and Σ is the population parameter covariance matrix. The core parameter of this layer is (μ,Σ), which is used to reflect the correlation between individual parameters and population parameters.

[0029] 3. Third layer (population parameter layer): Let the population parameter mean μ = Combining prior literature information and clinical data on individualized drug administration, The expression is:

[0030] ,

[0031] ;

[0032] in, The creatinine clearance rate reflects the patient's renal function and is a core individual characteristic in personalized medication. Age refers to the patient's individual age. Follows a mean of 0 and a variance of The normal distribution (k=1,2,3) The random effects term for the population parameter is used to reflect individual differences within the population, and is specifically expressed as: This is an individual difference correction term for the clearance rate CL population baseline value; Creatinine clearance rate Individual variability correction term for the effect of CL; This is the individual difference correction term for the distribution volume V, which represents the population baseline value.

[0033] 4. Fourth Layer (Prior Distribution Layer): Based on conjugate prior theory, combined with clinical experience and literature data in individualized dosing, the prior distribution of each parameter is set: the variance parameter of the random effect term is preset to...

[0034] =0.415、 =0.766、 =0.381,

[0035] The variance of the observation error term, σ, is 5.26 mg·L⁻¹ (adapted to the error range of clinical blood drug concentration observation); the inverse of the covariance matrix Σ follows a Wieshard distribution.

[0036]

[0037] The expected value is EΣ⁻¹=mΩ, meaning Σ follows an inverse Wissaud distribution, where m≥p=2 (p is the parameter dimension) to ensure Σ is invertible. Without clinical data support, a minimum information prior distribution is constructed, i.e., m=p=2 and Ω is a 2-order identity matrix I². When clinical observation data from n patients are introduced, then... The posterior distribution is as follows:

[0038] in:

[0039]

[0040]

[0041] in, Let S be the centered sample of this normal distribution, representing the deviation vector of individual parameters from the population mean, and let S be the sample covariance matrix, calculated from clinical data, expressed as follows:

[0042]

[0043] The formula for calculating the expected value is:

[0044]

[0045] Take the expected value of the posterior distribution as the covariance matrix. The estimated values ​​are used to ultimately estimate the population variance.

[0046] 5. Target Distribution Integration: Based on Bayesian theory, the above four-layer structure is integrated to construct the joint probability density function of the target probability distribution, i.e., the posterior distribution of the core parameters for individualized drug administration. Let the vector composed of numerical parameters be θ=( , If the target distribution is such that the posterior probability density function is:

[0047] Right now:

[0048]

[0049] because As a priori constant, we get:

[0050]

[0051] S12: Initialize the initial state of the Markov chain Based on prior information about the target distribution and the clinically normal ranges of parameters such as CL and V in individualized drug administration, an initial state is randomly generated within the parameter space. Record initial state information;

[0052] S13: Initialize sampling auxiliary parameters, including the acceptance rate reference threshold. (Values ​​range from 0.2 to 0.5, adapting to the acceptance rate requirements of individualized dosing parameter sampling), convergence judgment window size W (value range from 50 to 200), autocorrelation threshold. (Values ​​range from 0.1 to 0.3), and initialize the iteration counter t=0, and sample the sample set S= The chain state sequence X=[ ].

[0053] Furthermore, step S2 specifically includes:

[0054] S21: Based on the type and parameter dimensions of the target distribution in the individualized drug delivery scenario, select an appropriate fixed suggested distribution q( | It can adapt to the distribution characteristics of pharmacokinetic parameters;

[0055] S22: The parameters of the fixed suggested distribution are, based on existing literature, selected by this patent as a multivariate normal distribution, with the mean being the current chain state xt, and the covariance matrix Σ fixed as a preset identity matrix I or a preset constant matrix, similarly adapting to the fluctuation range of individualized dosing parameters, i.e., q( | )~N( Then, during a normal random walk, the generated candidate samples are:

[0056]

[0057] in For a single candidate sample, it refers to... The starting point, after taking one step according to the proposed distribution, is used as the reference point for the constraint boundary:

[0058]

[0059] The proposed distribution has a mean equal to the current state, satisfying the Markov property, and the covariance matrix is ​​chosen as a diagonal matrix for efficient computation.

[0060] S23: After determining the proposed distribution, keep the type and parameters of the proposed distribution unchanged throughout the entire sampling iteration process to meet the stability requirements of individualized dosing parameter sampling.

[0061] Furthermore, step S3 specifically includes:

[0062] In the traditional MCMC sampling process, a set of Euclidean distance-based constraints is introduced. With direction constraint set The dual-constraint mechanism maps candidate samples to [the target area] using a projection operator. ∩ This allows for the construction of constrained Markov chain state update paths.

[0063] S31: Iteration counter t=t+1, based on the current fixed proposal distribution q( | Choose the initial distribution and initial points. Initial sampling is performed to generate candidate samples. , Candidate sample vectors for core parameters of individualized drug administration are used as reference points for constraint boundaries;

[0064] S32: Define the parameter space as C (i.e., the reasonable value space of core parameters in individualized dosing, set based on the normal range of clinical parameters, such as the value range of CL being 0.1~10L / h, and the value range of V being 1~100L), based on the current sample, the previous sample, and the initial state. Construct two filter sets:

[0065] 321. Constructing a collection:

[0066] ={z∈C: }

[0067] Where ||·|| denotes the Euclidean norm of the vector, meaning that all vectors in the parameter space C are parallel to each other. The distance is less than or equal to the distance between The set of points z at which the distance is calculated is used to constrain the reasonableness of the distance of the samples and prevent the samples from deviating from the reasonable range of parameters;

[0068] 322. Constructing a set:

[0069] ={z∈C: ≥0}

[0070] Where <·,·> denotes the inner product of vectors, meaning that all vectors in the parameter space C that satisfy ( ) and vector ( The set of points z whose inner product is greater than or equal to 0 is used to constrain the iteration direction of the samples, thereby improving the convergence speed and consistency of the samples.

[0071] S323: Calculate the projection point:

[0072] ∩ ( )

[0073] in ∩ (·) indicates that the point will be... Projection to set With sets The intersection ( ∩ The projection operator on ) is found ∩ Zhongyu The closest point is used as the final sample after screening, through a fixed initial clinically reasonable point. This forces all candidate samples in all sampling steps to always converge toward a clinically reasonable region, thus preventing the Markov chain from diverging into a parameter space with no clinical significance during iteration.

[0074] like exist ∩ If the point is inside, then that point is taken as the projection result, that is... ;like Not here ∩ Inside, then calculate ∩ The closest point, the projection result is .

[0075] S324: Simultaneously, during the sampling process based on the projection operator, from the current state... Candidate samples were generated by proposing and projecting. The transition probability density is also a new proposed distribution. It can be represented as:

[0076]

[0077]

[0078]

[0079]

[0080] Among them, due to probability density normalization,

[0081] in Let be the Dirac function, representing a given and ,from Starting from there, you can only get one candidate.

[0082]

[0083] This density is a probabilistic representation of the deterministic mapping. Furthermore, since the projection operator is a deterministic mapping, its transition kernel still satisfies the detailed stationarity condition:

[0084]

[0085] in, For a complete transfer kernel including projection, For the target posterior distribution, and the transition kernel Represented as:

[0086]

[0087] Therefore, this method first performs projection preprocessing on the proposed samples. After preprocessing, the detailed balance condition is still satisfied, so the sampling can converge to the target distribution.

[0088] S33: After obtaining the new projection point, sampling is performed next. Since the parameters in this patent are numerical parameters, the resulting vector is θ=( , The covariance matrices Σ and Σ (2×2) have different parameter spaces. If they are in the same MCMC chain, it will lead to problems such as difficulty in designing the proposed distribution, slow convergence speed, and poor mixing effect. Therefore, this patent adopts a phased sampling method, first fixing the population covariance matrix Σ, and then using the MH algorithm to obtain the conditional posterior distribution. Adopt a six-dimensional numerical parameter vector θ=( , ); then, fixing the updated θ, based on the inverse Wieshard conjugate prior theory, from the conditional posterior distribution Sampling is performed to obtain Σ. The above sampling steps are repeated to obtain a sample chain. Converges on the joint posterior distribution of the target .

[0089] Let the vector of numerical parameters be θ=( , Based on the MH rule, the complete acceptance probability formula with proposal distribution is calculated.

[0090] α=min{ };

[0091] in, Let be the posterior probability density function of the target distribution. Substituting the conditional probability density function representing the proposal distribution into the newly derived proposal distribution above...

[0092]

[0093] The final acceptance probability can be obtained as follows:

[0094] α=min{ };

[0095] when and ,at this time Then the acceptance probability simplifies to α = min{ This indicates bidirectional transferability;

[0096] when At this time, in the molecule The function evaluates to 0, which necessarily means that the function will be rejected, indicating that... Cannot return via projection This indicates that the transfer is unidirectional, disrupting the detailed balance, therefore we reject the unidirectional transfer; or the denominator The function being 0 necessarily accepts the value, indicating that... Only by If we obtain the data, we set the acceptance probability to 1 to ensure that the chain can move forward and can completely collect the key areas.

[0097] S34: Update the current chain state, generate a uniformly distributed random number u in the interval [0,1]. If u < α, then accept. Update the state of the chain and add the new sample to the sample set. If u > α, then reject. The chain state remains unchanged.

[0098] S35: Fixed updated six-dimensional parameter vector Based on the inverse Wiesart conjugate prior theory, directly from

[0099]

[0100] Extracting a new covariance matrix Update the current chain state and add it to the covariance sample set.

[0101] S36: If To reject a candidate sample, return to the previous step and repeat steps S31 to S35 until the number of iterations t reaches the initial iteration threshold. .

[0102] S37: Pre-burning iteration. After collecting the sample set, samples that have not converged in the initial stage of the Markov chain need to be removed to avoid the influence of initial value deviation on the final effective samples.

[0103] First, set the base number of iterations during the pre-burn-in period. The value is taken as the initial iteration number threshold. 50%~80%; ensure sufficient iterations during the pre-burning period to allow the chain to escape the influence of the initial value. Execute iterations until t= The samples were screened during the pre-burning period (t=0 to t= All candidate samples and projected samples generated are discarded and do not participate in subsequent convergence judgment, autocorrelation removal, and statistical feature calculation. Only t> The subsequent chain state sequence is used as the effective convergence judgment sequence and enters the convergence judgment stage in step S4.

[0104] Furthermore, step S4 specifically includes:

[0105] S41: It adopts traditional convergence judgment indicators, including two core indicators, the potential scale reduction factor (PSRF) and the autocorrelation coefficient, to meet the convergence requirements of individualized dosing parameter sampling.

[0106] 1) Potential Scale Reduction Factor (PSRF): Calculate the PSRF value of the chain state sequence X. If PSRF ≤ 1.1, the chain is considered to be stationary in this dimension, that is, the parameter samples are stable.

[0107] 2) Autocorrelation coefficient: Calculate the autocorrelation coefficient of each dimension parameter in the chain state sequence X. If the autocorrelation coefficient of all dimensions is ≤ when the lag order is W... If the mixing property of the chain is satisfactory, the independence of the samples is considered to be good.

[0108] S42: If both core metrics meet the requirements, the chain is considered converged, and proceed to step S5; if at least one metric does not meet the requirements, the number of iterations is increased (each iteration increases by 100%). / 10), return to step S3 to continue iterating until the convergence condition is met; if the number of iterations reaches the maximum iteration threshold. ( =10 If convergence is still not achieved, a convergence warning will be output and the process will proceed to step S5 (to avoid sampling failure due to convergence issues, which would affect the efficiency of individualized dosing design).

[0109] Furthermore, step S5 specifically includes:

[0110] S51: Autocorrelation removal, determining the sampling step size k based on the autocorrelation coefficient (k is the step size for reducing the autocorrelation coefficient to...). The following minimum lag order is used to draw one sample from the sample set S every k samples to obtain the effective sample set after deautocorrelation. (i.e., a valid sample of core parameters for individualized drug administration).

[0111] S52: If the valid sample set The number of samples is less than the preset minimum number of valid samples. ( =0.5 If N), then return to step S3, increase the number of iterations, and supplement the sampled samples until N). The sample size meets the requirements (ensuring that the sample size is sufficient to support the estimation of individualized dosing parameters).

[0112] Furthermore, step S6 specifically includes:

[0113] S61: Calculate the effective sample set The statistical characteristics, including the mean, variance, median, and quantiles of pharmacokinetic parameters in each dimension, as well as the joint probability distribution histogram of the samples, provide data support for the optimization of individualized dosing regimens.

[0114] S62: Sample validity verification. The consistency between the sample and the target distribution is verified using the marginal distribution KS test and the joint structure KDE fit test. KS tests are performed on each marginal distribution of the joint posterior. If the p-value of all marginal distributions is greater than 0.05, the marginal structure verification is successful. Next, a subset of core clinical parameters is selected, and kernel density estimation (KDE) is used to fit the joint probability density of the sampled sample and the theoretical joint posterior simulated sample, respectively, calculating their KL divergence. If the KL divergence is less than 0.1, the joint structure verification is successful. If both steps pass, the valid sample set is determined to be consistent with the target joint posterior distribution, and the verification is valid; otherwise, a verification warning is output, and the process returns to step S3 for resampling.

[0115] S63: Output the set of valid samples The statistical characteristics and validation results provide core parameter support for the design of individualized dosing regimens and complete the entire sampling process.

[0116] Beneficial effects

[0117] The sampling method and system for personalized drug delivery provided by this invention have the following advantages compared to existing technologies:

[0118] 1. Strong scenario adaptability: With personalized drug administration as the core background, the hierarchical analysis method of ppk+MAPB algorithm is used to construct a multi-level nested target distribution, which accurately adapts to the hierarchical relationship of "observation data-individual parameters-population parameters-prior information". This solves the pain point that traditional MCMC sampling cannot adapt to the complex parameter distribution of personalized drug administration. It can be directly applied to personalized drug administration related scenarios such as pharmacokinetic parameter estimation and blood drug concentration prediction.

[0119] 2. Significantly improved sample quality: The addition of dual-set screening and projection operations in step S3 further purifies the samples through distance constraints and inner product constraints, filtering out invalid samples that deviate from the reasonable range of parameters or have excessively high autocorrelation, making the final samples more consistent with the distribution of individualized dosing targets, improving the accuracy of parameter estimation, and providing reliable data support for the optimization of individualized dosing regimens.

[0120] 3. Eliminate abnormal parameters without clinical significance: By setting a clinically reasonable range for the parameter space C, combined with a set of distance constraints. By filtering out candidate samples that deviate from a reasonable range, abnormal parameters that meet the probability distribution but do not conform to physiological logic are directly removed from the sampling stage. This solves the problem that traditional MCMC sampling results are often statistically valid but clinically invalid, and avoids calculation errors in subsequent dosing regimens due to abnormal parameters.

[0121] 4. Iterative directional constraints to achieve directional contraction of the sampling space: through a set of directional constraints The inner product constraint forces the Markov chain to always converge toward the clinically reasonable region, solving the problem of repeated oscillations in traditional random walks, making the sampling iteration more directional and effectively reducing the number of invalid iterations.

[0122] 5. The projection operator ensures that samples are always within the clinically reasonable domain: the projection operator maps candidate samples to... Even if there are slight deviations in the initial values, the samples can be pulled back to the clinically reasonable parameter space through projection, while avoiding the divergence of the Markov chain with iteration, thus greatly improving the overall quality of the sampled samples.

[0123] 6. Reduce invalid sampling ratio and improve computational efficiency: Dual-constraint screening filters out invalid candidate samples that are likely to be rejected in advance, reducing the repetitive "generate-reject-regenerate" operation in traditional sampling, significantly reducing computational resource consumption and shortening sampling time.

[0124] 7. Solving the sampling adaptation problem of different parameter spaces: To address the parameter space difference between the numerical parameter vector θ and the covariance matrix Σ, an alternating strategy of "first fixing Σ and sampling θ, then fixing θ and sampling Σ" is adopted. This avoids the problems of difficult distribution design and poor mixing effect in traditional single-chain sampling, allowing the sampling of both types of parameters to adapt to their own distribution characteristics and improving the stability of sampling. Attached Figure Description Figure 1 This is a flowchart of the method of this patent. Detailed Implementation

[0127] The present invention will be further described in detail below with reference to specific embodiments. This embodiment is based on the individualized dosing scenario that conforms to the first-order elimination and one-compartment model under the background of vancomycin intravenous infusion administration. It is based on population pharmacokinetics and the maximum a posteriori Bayesian estimation algorithm, and assumes that relevant pharmacokinetic parameters are sampled, which reflects the core improvement of the present invention and adapts to the clinical needs of individualized dosing.

[0128] Example of a hierarchical mathematical model for the ppk+MAPB algorithm

[0129] S1: Determine the prior data conditions and the target posterior distribution:

[0130] ①Based on cutting-edge literature, statistical clinical data, and the basic assumptions of models such as PPK, we conducted a study on n actual patients to obtain individual... Actual measured data such as age, weight, and renal function CLcr were used to calculate the sample covariance matrix, the mean and variance of CL and V for the population and the sample.

[0131] ②Based on pharmacokinetic principles, the blood drug concentration after continuous administration The formula is

[0132]

[0133] Used to calculate predicted blood drug concentrations for individuals;

[0134] ③Let obj is the observed value of individual i at time j, specifically the observed blood drug concentration in a personalized drug administration scenario. `pre` is the model-predicted blood drug concentration at time `j` for individual `i`, derived from ②. We obtain the result based on individual random effects. = The relationship between observed and predicted values ​​is obtained:

[0135] ln obj=ln pre+ ,

[0136] in It is a random error term and:

[0137]

[0138] According to the literature Therefore, we get:

[0139] ln obj

[0140] ④ Based on the fundamental assumptions of the PPK model, if multiple pharmacokinetic parameters in the population follow a multidimensional normal distribution, then:

[0141]

[0142] Let the population parameter mean μ = ( According to the ppk covariate model, there is a quantitative relationship between the population pharmacokinetic parameter mean and the individual covariates, namely:

[0143]

[0144] Based on existing baseline drug values ​​and the relationships between various factors, and considering clinical trials and literature, we hypothesize that:

[0145] ,

[0146] ;

[0147] At the same time, set

[0148]

[0149] The variance parameter of the random effects term is set as =0.415、 =0.766、 =0.381,

[0150] The covariance matrix is ​​in the form of Specifically, it is obtained from the Bayesian estimation in ⑤.

[0151] ⑤ Based on the conjugate prior theory, combined with clinical experience and literature data on individualized drug administration, ④ the inverse matrix of the covariance matrix Σ follows a Wieshard distribution:

[0152]

[0153] expect =mΩ, meaning Σ follows an inverse Wissaud distribution, where m ≥ p = 2 (p is the parameter dimension) to ensure Σ is invertible. Without clinical data support, a minimum information prior distribution is constructed, i.e., m = p = 2 and Ω is a 2-order identity matrix. When clinical observation data from n patients are introduced, the following is obtained: The posterior distribution is as follows:

[0154]

[0155] in:

[0156]

[0157]

[0158] in, Let S be the centered sample of this normal distribution, representing the deviation vector of individual parameters from the population mean, and let S be the sample covariance matrix, calculated from clinical data, expressed as follows:

[0159]

[0160] The formula for calculating the expected value is:

[0161]

[0162] Take the expected value of the posterior distribution as the covariance matrix. The estimated values ​​are used to ultimately estimate the population variance.

[0163] ⑥ Target Distribution Integration: Based on Bayesian theory, the above four-layer structure is integrated to construct the joint probability density function of the target probability distribution, i.e., the posterior distribution of the core parameters for individualized drug administration. Let the vector composed of numerical parameters be θ=( , If the target distribution is such that the posterior probability density function is:

[0164] Right now:

[0165]

[0166] because As a priori constant, we get:

[0167]

[0168] Based on the probability density functions and the fact that only the relative density of the posterior distribution is needed, the parameter-independent constant term can be discarded for simplified calculation, ultimately yielding:

[0169]

[0170] in:

[0171]

[0172]

[0173]

[0174]

[0175] S2: Initialize the initial state of the Markov chain. Based on the prior information of the target distribution and the clinically normal ranges of parameters such as CL and V in individualized drug administration, define the initial distribution in the parameter space, and randomly sample a point from the initial distribution. As the initial state of the chain, record the initial state information;

[0176] S3: Initialize sampling auxiliary parameters, including the acceptance rate reference threshold. (Values ​​range from 0.2 to 0.5, adapting to the acceptance rate requirements of individualized dosing parameter sampling), convergence judgment window size W (value range from 50 to 200), autocorrelation threshold. (Values ​​range from 0.1 to 0.3), and initialize the iteration counter t=0, and sample the sample set S= The chain state sequence X=[ ].

[0177] Specific implementation steps of the sampling section

[0178] S1: Sampling initialization

[0179] S11: Determine the core sampling parameters. The target distribution is a multi-level nested continuous distribution constructed based on Bayesian hierarchical analysis, with parameter dimension d=6, and θ=( , The target probability density function is the joint posterior distribution described above, the total number of samples is N=1000, and the initial iteration number threshold is... =1000, maximum iteration threshold =10000;

[0180] S12: Randomly initialize the initial state Within the parameter space C, the initial state is randomly generated. =(2.0,25.0,0.1,0.2,0.1,5.26) (CL=2.0L / h, V=25.0L, =0.1, =0.2, =0.1, =5.26), which is within the reasonable range of clinical parameters;

[0181] S13: Initialize auxiliary parameters, acceptance rate reference threshold =0.3, convergence judgment window size W=100, autocorrelation threshold =0.2, iteration counter t=0, sample set S= The chain state sequence X = [(2.0, 25.0, 0.1, 0.2, 0.1, 5.26)].

[0182] S2: Determine the fixed recommended distribution

[0183] S21: Select a multivariate normal distribution as the proposed distribution q( | )~N( Let I be a 6×6 identity matrix (adapted to 6-dimensional parameter sampling). Then, during a normal random walk, the generated candidate samples are... , as the reference point for the constraint boundary:

[0184] ;

[0185] S22: At this point, all the required prerequisites are obtained: the parameter vector is θ=( , The joint posterior distribution of the target and the initial state of the Markov chain generated based on the parameter space C. Proposal distribution, auxiliary parameters (acceptance rate) =0.3, convergence judgment window size W=100, autocorrelation threshold =0.2, iteration counter n=0, sample set S= The chain state sequence X = [(2.0, 25.0, 0.1, 0.2, 0.1, I2)] and the parameter space C range and projection operator.

[0186] S3: Markov chain iterative sampling

[0187] S31: n iterates starting from 1, and each iteration is based on the current fixed proposal distribution q( | ), generate candidate samples That is, candidate values ​​of parameters such as CL, V, η1, η2, η3, etc., ensuring that the candidate values ​​are within the parameter space C, and serving as reference points for the constraint boundaries;

[0188] S32: Order For the newly generated constraint boundary reference points, perform the following additional filtering step; after filtering, if u > α, reject. Repeat step S33 until an acceptable candidate sample is obtained; if u ≤ α, accept. Then, we begin the next iteration.

[0189] S321: Define the parameter space as C (CL∈[0.5~5L / h], V∈[10~50L]), based on the current sample, the previous sample, and the initial state. Construct two filter sets:

[0190] 1. Set ={z∈C: }, where ||·|| represents the Euclidean norm of the vector, i.e., all vectors in the parameter space C that are parallel to each other. The distance is less than or equal to the distance between The set of points z at which distance is given; where:

[0191]

[0192] in, Let be the k-th element of the vector. This formula is essentially about using in 6-dimensional space... and The perpendicular bisector of the connecting line forms the boundary of the hyperplane, close to... The half-space.

[0193] S322: Set ={z∈C: ≥0}, where <·,·> denote the inner product of vectors, that is, all vectors in the parameter space C that satisfy ( ) and vector ( The set of points z whose inner product is greater than or equal to 0;

[0194]

[0195] The formula is essentially a pass / fail. And with vectors ( A vertical hyperplane with a non-negative inner product of half-spaces.

[0196] S323: Calculation ∩ ( ),in ∩ (·) indicates the initial state Projection to set With sets The projection operator on the intersection, i.e., finding ∩ Zhongyu The closest point is selected as the final filtered sample; that is:

[0197]

[0198] The formula for the squared Euclidean distance in six dimensions is:

[0199]

[0200] like If it is not within the projection area, then project it to... Solve , and then Projecting the initial point to Solving for the final result If it is within the projection area, then take it directly. As

[0201] S324: Update the current chain state ,like If a candidate sample is accepted, then... Add the sample set S, and then according to... Perform a normal random walk to generate the next step. ;like To reject candidate samples, then not Add to the numerical parameter sample set Simultaneously take the previous step Repeat the above steps with a normal random walk until the iteration ends.

[0202] S325: Simultaneously, during the sampling process based on the projection operator, from the current state... Candidate samples were generated by proposing and projecting. transition probability density It can be represented as:

[0203]

[0204]

[0205]

[0206]

[0207] Among them, due to probability density normalization,

[0208] in Let be the Dirac function, representing a given and Projection results It is uniquely determined that this density is a probabilistic representation of a deterministic mapping. Furthermore, since the projection operator is a deterministic mapping, its transition kernel still satisfies the detailed stationarity condition:

[0209]

[0210] in, For a complete transfer kernel including projection, For the target posterior distribution, and the transition kernel Represented as:

[0211]

[0212] Therefore, this method first performs projection preprocessing on the proposed samples. After preprocessing, the detailed balance condition is still satisfied, so the sampling can converge to the target distribution.

[0213] S33: Calculate the probability of acceptance

[0214] α=min{ };

[0215] Generate uniformly distributed random numbers u in the interval [0,1] and accept them according to probability to satisfy the detailed balance condition.

[0216] S34: Accept newly generated samples After that, it was obtained Fixed updated six-dimensional parameter vector Based on the inverse Wiesart conjugate prior theory, directly from

[0217]

[0218] Extracting a new covariance matrix Update the current chain state and add it to the covariance sample set. ;

[0219] S35: Will Add X, iterate repeatedly until t=1000, and then proceed to the convergence judgment stage.

[0220] S36: Pre-burning iteration. After collecting all samples, samples that have not converged in the initial stage of the Markov chain need to be removed to avoid the impact of initial value bias on the final valid samples. First, set the basic number of pre-burning iterations. The value is taken as the initial iteration number threshold. 50%~80%, in this embodiment, we take =0.6× ,because =1000, then =600, ensuring that there are enough iterations during the pre-burning period to allow the chain to get rid of the influence of the initial value.

[0221] Execute iterations up to t= The samples were screened during the pre-burning period (t=0 to t= All candidate samples and projected samples generated are discarded and do not participate in subsequent convergence judgment, autocorrelation removal, and statistical feature calculation. Only t> The subsequent chain state sequence is used as the effective convergence judgment sequence and enters the convergence judgment stage in step S4.

[0222] S4: Convergence test

[0223] S41: Calculate the PSRF value and autocorrelation coefficient of each dimension of the chain state sequence X. If PSRF ≤ 1.1 and autocorrelation coefficient ≤ 0.2 when lag order is 100, then convergence is determined.

[0224] In this embodiment, when t=1000, PSRF=1.08, and the autocorrelation coefficients of all dimensions are ≤0.18, which meets the convergence condition, so proceed to step S5.

[0225] S5: Sample Screening

[0226] S51: Calculate the minimum lag order k=10 for the autocorrelation coefficient to drop below 0.2, and draw one sample every 10 samples to obtain the effective sample set. Its expression is:

[0227]

[0228] Where M is the number of valid samples, Let be the parameter vector of the kth valid sample, and let θ be the parameter vector in the target distribution. , ),Right now =( , ).in, Let be the clearance rate of individual i in the k-th valid sample; The volume of distribution for individual i in the k-th valid sample (unit: L); This represents the random effects term of the population parameters in the k-th valid sample; The variance of the observation error term is σ = 5.26 mg·L⁻¹. In this embodiment, the sample size is 610, which satisfies... =500 requirement;

[0229] S6: Sampling Result Output and Verification

[0230] S61: Calculation The statistical characteristics are used to calculate the mean, which reflects the population average level of the parameter and is used as the basis for setting individualized dosing dosages:

[0231] ,

[0232]

[0233] in The effective sample mean of CL. V is the effective sample mean, and the effective sample size M = 610; then the variance is calculated to reflect individual differences in the parameters, which is used for individualized adjustment of individualized dosing regimens:

[0234] ,

[0235]

[0236] in, Let CL be the effective sample variance. Let V be the effective sample variance; then calculate the statistical characteristics of the random effects term and the covariance matrix:

[0237] ,

[0238]

[0239] in, For random effects The effective sample mean, Covariance matrix The effective sample mean is used to analyze the correlation between individual parameters within a population.

[0240] In this embodiment, the mean CL was 2.1 L / h with a variance of 0.32, and the mean V was 24.8 L with a variance of 4.5. The mean values ​​are all close to 0, which conforms to the statistical characteristics of the target distribution and is suitable for the distribution pattern of individualized dosing parameters.

[0241] S62: Sample validity verification. The consistency between the sample and the target distribution is verified using the marginal distribution KS test and the joint structure KDE fit test. KS tests are performed on each marginal distribution of the joint posterior. If the p-value of all marginal distributions is greater than 0.05, the marginal structure verification is successful. Next, a subset of core clinical parameters is selected, and kernel density estimation (KDE) is used to fit the joint probability density of the sampled sample and the theoretical joint posterior simulated sample, respectively, calculating their KL divergence. If the KL divergence is less than 0.1, the joint structure verification is successful. If both steps pass, the valid sample set is determined to be consistent with the target joint posterior distribution, and the verification is valid; otherwise, a verification warning is output, and the process returns to step S3 for resampling.

[0242] S63: Output Statistical characteristics and validation results, based on effective samples of CL and V, combined with patient age, Based on individual characteristics, optimize the dosage (e.g., adjust the dosage to 500 mg, with a dosing interval of 12 h) and complete the sampling.

[0243] S64: The following data results were obtained through simulation:

[0244] Table 1

[0245]

[0246]

Claims

1. A Markov chain Monte Carlo sampling method based on a double-constraint set projection mechanism, suitable for personalized drug delivery scenarios, characterized in that, Includes the following steps: S1: Sampling initialization, determine the core sampling parameters in the individualized dosing scenario, and use the population pharmacokinetic and maximum a posteriori Bayesian estimation hierarchical analysis (ppk+MAPB) algorithm to construct a four-layer nested target probability distribution model of observation layer - individual parameter layer - population parameter layer - prior distribution layer, initialize the initial state of the Markov chain, and initialize the sampling auxiliary parameters at the same time; S2: Determine the proposal distribution. Based on the distribution characteristics of individualized dosing parameters, select a multivariate normal distribution as the proposal distribution. Set the mean of the proposal distribution to the current chain state and the covariance matrix to be a diagonal matrix. S3: Markov chain iterative sampling, based on the Metropolis-Hastings criterion, combined with the proposal distribution to generate candidate samples, and introducing a set of Euclidean distance constraints. With direction constraint set The dual-constraint mechanism maps candidate samples to [the target area] using a projection operator. ∩ After obtaining the filtered samples, the covariance matrix Σ is first fixed, and then the numerical parameter vectors are sampled. Then fix The phased sampling strategy of sampling Σ combines the acceptance probability formula to determine whether to accept new samples, update the chain state and remove samples that have not converged during the pre-burning period; S4: Convergence judgment, using the latent scale reduction factor and autocorrelation coefficient as core indicators to monitor the stability of the chain. If the indicators do not meet the requirements, the number of iterations is increased and resampling is performed until the convergence condition is met or the maximum iteration threshold is reached. S5: Sample screening. The sampling step size is determined based on the autocorrelation coefficient of the converged sample, and autocorrelation removal is performed to obtain effective samples. If the number of effective samples is insufficient, additional sampling is performed. S6: Sampling results output and verification. Calculate the statistical characteristics of the effective samples, and use the marginal distribution KS test and the joint structure KDE goodness test to verify the consistency between the samples and the target distribution. After the verification is passed, output the effective samples and statistical characteristics to provide data support for the optimization of individualized dosing regimens.

2. The sampling method according to claim 1, characterized in that, The four-layer nested target probability distribution model described in step S1 is constructed as follows: Observation layer: Establish the logarithmic linear relationship between observed blood drug concentrations and model predictions, ln obj=ln pre+ , ln pre is determined by the individual's core pharmacokinetic parameter clearance rate Distributed volume Determined by a preset concentration calculation function h(·); defined as h(·) = , Individual parameter layer: Defines the individual's core pharmacokinetic parameter vector. , The population pharmacokinetic parameters are the overall mean. The group parameter covariance matrix; Population parameter layer: combined with patient creatinine clearance rate Individual characteristics such as age are used to construct the mean of group parameters. = The quantitative calculation formula is given, and the random effects term of the population parameter is set. ; Prior distribution layer: Based on conjugate prior theory, a priori distribution layer is defined. ,Right now Following an inverse Wissaud distribution, clinical observation data were introduced to obtain... The posterior distribution is given, and the variance of the random effects term is also specified. Prior constants for clinical adaptation; ultimately, based on Bayesian theory, a four-layer structure is integrated to construct the joint probability density function of the target probability distribution, with a numerical parameter vector θ=( , ).

3. The sampling method according to claim 1, characterized in that, The sampling auxiliary parameters mentioned in step S1 include: acceptance rate reference threshold. The value ranges from 0.2 to 0.5; the convergence judgment window size W ranges from 50 to 200; the autocorrelation threshold... The value ranges from 0.1 to 0.3; at the same time, the iteration counter, the sample set, and the chain state sequence are initialized.

4. The sampling method according to claim 1, characterized in that, The Euclidean distance constraint set mentioned in step S3 and direction constraint set The definition of is: ={z∈C:|| |≤|| ||}, where C represents the clinically reasonable value space for the core parameters of individualized drug administration. As candidate samples, Let ||·|| represent the current chain state, and ||·|| denote the Euclidean norm of the vector. ={z∈C:< , >≥0}, Let P be the initial state of the Markov chain, and <·,·> denote the inner product of vectors; the projection operator P is... ∩ (·) Initial state Mapped to ∩ Zhongyu The nearest point yields the filtered samples. .

5. The sampling method according to claim 4, characterized in that, In the clinically reasonable value space C of the core parameters mentioned in step S3, the clearance rate CL ranges from 0.1 to 10 L / h, and the distribution volume V ranges from 1 to 100 L.

6. The sampling method according to claim 1, characterized in that, The specific operation of the phased sampling strategy described in step S3 is as follows: With a fixed covariance matrix Σ, the conditional posterior distribution is obtained from the MH algorithm. Sampling a six-dimensional numerical parameter vector θ=( ,η1,η2,η3, ); With the updated θ fixed, based on the inverse Wieshard conjugate prior theory, from the conditional posterior distribution... Sample covariance matrix Σ; Alternately perform the above two steps to make the sample chain... Converges on the joint posterior distribution of the target .

7. The sampling method according to claim 1, characterized in that, The acceptance probability formula mentioned in step S3 is: α = min{ },in Let be the posterior probability density function of the target distribution. For the Dirac function, The projection operator generates a uniform random number u in the interval [0,1]. If u < α, then the input is accepted. Otherwise, reject the chain and keep the chain state unchanged.

8. The sampling method according to claim 1, characterized in that, The method for eliminating samples during the pre-burning period in step S3 is as follows: set the basic number of iterations for the pre-burning period. Initial iteration number threshold 50%~80%, excluding t=0 to t= All samples generated within the timeframe, only those with t> The subsequent chain state sequence is used as the effective convergence judgment sequence.

9. The sampling method according to claim 1, characterized in that, The convergence criterion mentioned in step S4 is: the potential scale reduction factor PSRF ≤ 1.1, and the autocorrelation coefficients of all dimensional parameters are ≤ 1.1 when the lag order is W. ; If the number of iterations reaches the maximum iteration threshold =10× If convergence is still not achieved, a convergence warning will be output and the sample selection step will be initiated.

10. The sampling method according to claim 1, characterized in that, The standard for validating the sample validity in step S6 is: perform KS test on each marginal distribution of the joint posterior, and if the p-value of the KS test for all marginal distributions is greater than 0.05, the marginal structure verification is passed; Select a subset of core clinical parameters and use kernel density estimation (KDE) to fit the joint probability density of the sampled samples and the theoretical joint posterior simulated samples. If the KL divergence between the two is less than 0.1, the joint structure validation is passed. If both steps are passed, the sample is determined to be consistent with the target joint posterior distribution. Otherwise, a validation warning is output and resampling is performed.

Images