A complex long-acting formulation steady-state pharmacokinetics extrapolation equivalence prediction method

Abstract

The application discloses a complex long-acting preparation steady-state pharmacokinetics extrapolation equivalence prediction method, relates to the steady-state bioequivalence technical field, and comprises the following steps: collecting and unifying concentration-time, a drug administration scheme and a covariate, generating a standardized time, an event marker code and a covariate dictionary, setting a deletion weight and a population alignment weight; setting a candidate structure based on population pharmacokinetics, estimating a model and a covariate effect through a quality threshold determination with a weight; constructing a virtual population and performing multiple drug administration simulation to determine the reaching of stability according to the relative change of consecutive two trough concentrations, generating a peak-trough sampling time window; performing non-compartment analysis in the steady-state interval, calculating a dosing interval area, a steady-state peak value and a steady-state trough value, calculating a geometric mean ratio and a confidence interval of two preparations, carrying out efficacy-sample size linkage, and outputting a recommended design; the method can reduce the test burden, shorten the cycle and improve the decision transparency.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of steady-state bioequivalence, in particular to a steady-state pharmacokinetic extrapolation method for predicting the equivalence of complex long-acting preparations. BACKGROUND

[0002] As an important type of high-end pharmaceutical preparations, complex long-acting preparations are commonly used in long-acting injection, sustained-release and implantation drug delivery systems, and have the advantages of long-lasting effect, reduced frequency of administration, and improved compliance. However, the clinical conversion rate and success rate of such preparations have been low for a long time. The current steady-state pharmacokinetic and bioequivalence studies still mainly rely on the traditional clinical protocol of multiple dosing to reach steady state and intensive blood sampling before and after reaching steady state. Due to the long half-life, slow absorption, and coexistence of inter-individual and intra-individual variability of such preparations, the steady-state period is long. The cross-over design is affected by compliance and intra-individual variability, and the parallel design is more sensitive to sample size and dropout. The test time and cost increase dramatically, and the cost of a single study is at least two million yuan. Moreover, any minor process change often requires repeated testing, and the risk of failure is difficult to predict in a timely manner.

[0003] From a technical point of view, it is urgently needed in clinical practice to evaluate the area under the curve, steady-state peak concentration and steady-state trough concentration within the dosing interval under the background of long-term medication, in order to support the prediction of efficacy and safety and individualized dosing.

[0004] However, the traditional approach mainly relies on non-compartmental analysis to calculate the observed data after reaching steady state, which cannot reliably extrapolate between different dosing scenarios, and it is difficult to scientifically determine the number of doses required to reach steady state and the key sampling window. The effects of preparation properties, physiological and pathological differences, and genetic differences are difficult to quantify, the utilization efficiency of existing data is low, and repeated or expanded trials have become the norm, increasing the ethical burden.

[0005] Although in theory and practice, quantitative pharmacology and population pharmacokinetic modeling have been applied in many aspects of drug development, there is still a lack of systematic methods and operational specifications widely adopted in the industry in the context of steady-state prediction and equivalence evaluation of complex long-acting preparations, making it difficult for researchers to optimize the protocol, control the sample size, and manage the risks in a timely and accurate manner.

[0006] Therefore, the current problem is: in the research and development and equivalence evaluation scene of complex long-acting preparations, under the influence of factors such as long half-life, slow reaching steady state, high variability, strong test burden and ethical constraints, the existing analysis paradigm based on observation after reaching steady state cannot reliably infer the steady-state exposure under long-term administration from the information before reaching steady state under the condition of only single or few doses and small sample clinical data, further cannot scientifically determine the number of doses required to reach steady state and the steady-state sampling time window, and it is also difficult to make a prospective prediction of the equivalence of two preparations at the steady-state level and quantify the failure risk, so that the research design, resource allocation and process change decision lack verifiable technical basis. SUMMARY

[0007] (I) Technical problems solved

[0008] In view of the deficiencies of the prior art, the present application provides a complex long-acting preparation steady-state pharmacokinetic extrapolation equivalence prediction method, which comprises setting a candidate structure based on population pharmacokinetics, estimating a model determined by a quality threshold and a covariate effect with weight; construct a virtual population and perform multiple dose simulation to determine steady state by the relative change of consecutive two trough concentrations, generate peak-trough sampling time window; non-compartment analysis in steady-state interval, calculate dosing interval area, steady-state peak and steady-state trough, calculate geometric mean ratio and confidence interval of two preparations and carry out power-sample size linkage, output recommended design; the method can reduce the test burden, shorten the period and improve the decision transparency; solve the technical problems recorded in the background art.

[0009] (II) Technical solutions

[0010] In order to achieve the above purpose, the present application is realized by the following technical solutions:

[0011] A complex long-acting preparation steady-state pharmacokinetic extrapolation equivalence prediction method, comprising: collecting single / few dose human concentration-time, dosing regimen and covariates and unifying the caliber, generating a deletion weight for low limit observation mark and calculating a population alignment weight, outputting modeling input data set, data processing record and reference model item;

[0012] Based on the modeling input data set, set the candidate structure containing the number of compartments, absorption and elimination form, estimate the parameters by using nonlinear mixed effects combined with two types of weights, determine the final model and covariate effect according to the quality threshold;

[0013] Based on the final model and covariate effects, a virtual population is generated according to the population covariate distribution, and multiple dosing simulations are performed. The threshold of relative change in trough concentration before two adjacent dosings is used to determine the steady state, and the peak and trough sampling time windows are determined according to the rate of change threshold within the steady state interval. Within the complete dosing interval after the steady state is reached, the dosing interval area, steady state peak value, and steady state trough value are calculated. The geometric mean ratio of the steady state indices of the subjects and the reference is calculated and determined according to the confidence and equivalence intervals. The efficacy is simulated under the sample size and enumeration of inter-individual differences, and the recommended design is output.

[0014] Furthermore, a covariate dictionary is established and units and codes are standardized. Population alignment weights are determined by the ratio of covariate marginal density between the target population and the sample, and interval pruning is performed. Records below the quantitative lower limit are labeled with censoring weights. A modeling input dataset with version fingerprints and data processing records are generated.

[0015] Furthermore, the candidate structure includes parallel zero-order release and first-order absorption with hysteresis time, and inter-chamber exchange between the central chamber and the peripheral chamber; parameter estimation is performed using weighted likelihood with incorporation population alignment weights and censoring weights; observation errors are modeled on a logarithmic scale; and the final model and covariate effect function are output.

[0016] Furthermore, the quality threshold includes at least the following: generating a goodness-of-fit plot, a time-stratified prediction validation plot, and bootstrap parameter reestimation; performing covariate screening according to forward and backward thresholds and fixing the random seed; setting an upper limit on the uncertainty of the main parameters and recording it to form an evidence package for passing the quality threshold.

[0017] Furthermore, the covariate correlation structure is reconstructed based on the covariate dictionary, and a list of virtual population parameters is generated; an event-time series is generated according to the dosing interval and start time and aligned with the standardized time and event identifier code; and the random seed and event list are fixed by document logging.

[0018] Furthermore, the determination of stabilization is based on the relative change in trough concentration before two consecutive administrations being lower than a preset threshold; within the complete dosing interval after stabilization, the sampling time window for peak and trough values ​​is determined based on the proportion of local change rate of the concentration curve not exceeding the reference slope, and a list of sampling windows with timestamps is generated.

[0019] Furthermore, within the first complete dosing interval after stabilization, a sampling grid is constructed with a sampling window list as the priority constraint. The area within the dosing interval is calculated using trapezoidal integrals, and steady-state peak values ​​and steady-state valley values ​​are extracted. Peak values ​​and valley values ​​are searched and determined one by one within the corresponding time window, and the corresponding time information is recorded.

[0020] Furthermore, the geometric mean ratio and confidence interval of the steady-state indices of the test formulation and the reference formulation are calculated on a logarithmic scale, and equivalence is determined by the inclusion of the preset equivalence interval; wherein, the degrees of freedom and design type are matched and fixed in the equivalence determination report.

[0021] Furthermore, efficacy simulations were conducted on a two-dimensional enumeration grid of inter-individual differences and sample size, with the number of repetitions and random seeds fixed. Among combinations that meet the efficacy threshold, engineering constraints of the total number of blood collections and the study period were superimposed to screen the study design, and the recommended design was output in document form.

[0022] Furthermore, the entire process uses the modeling input dataset, reference model item list, data processing records, stability determination report, sampling window list, steady-state index list, and equivalence determination report as the sole input and output criteria, and maintains a unique mapping between terms and variable identifiers, all of which are written into version fingerprints and timestamps.

[0023] (III) Beneficial Effects

[0024] This invention provides a method for predicting the equivalence of steady-state pharmacokinetic extrapolation in complex long-acting formulations, which has the following beneficial effects:

[0025] Starting with data access and availability preparation, the concentration-time, dosing regimen, and covariate definitions are standardized, and standardized time, event flag codes, and covariate dictionaries are generated. Censoring weights are set for lower limit observations, and population alignment weights are set for population differences. Data processing records are solidified to provide an entry point for subsequent modeling and extrapolation. In population pharmacokinetic modeling and quality threshold validation, candidate structures cover zero-order and first-order absorption and include hysteresis and central and peripheral compartment exchange. Weighted likelihood incorporates the two types of weights into the estimation, and covariates are entered into the parameters in the form of functions. The model is solidified through goodness-of-fit, prediction verification, and bootstrap thresholding to ensure consistency of definitions.

[0026] In virtual population simulation, the relevant structure is reconstructed based on the covariate dictionary, and a list of virtual population parameters is generated. Multiple dosing simulations are carried out by aligning event-time series with standardized time and event flag codes. The steady-state threshold is determined by using the relative change in trough concentration before two consecutive dosings to determine the number of dosings and the steady-state starting point. Within the same steady-state interval, sampling time windows are automatically generated around the peak and trough values ​​using the rate of change threshold, and timestamps are bound to file traces. The sampling location and time tolerance range are unified to reduce the inconsistency of peak-trough offset and ensure that subsequent calculations and judgments are performed on the same time baseline.

[0027] In the non-compartmental analysis of the steady-state interval, the dosing interval area is calculated according to a uniform sampling grid, and steady-state peak and valley values ​​are extracted. The geometric mean ratio and confidence interval are calculated for the steady-state indices of the test and reference, so that the equivalence determination is based solely on steady-state information, avoiding interference from the transition segment and forming a decision chain. The steady-state equivalence determination is linked with efficacy-sample size enumeration, with fixed degrees of freedom and random seeds, and engineering constraints of total blood sampling and study period are superimposed to output the recommended study design. The entire process uses version fingerprints and timestamps to solidify input-processing-output, facilitating auditing and reuse. Attached Figure Description

[0028] Figure 1 This is a schematic diagram of the steady-state pharmacokinetic extrapolation equivalence prediction method for complex long-acting formulations according to the present invention. Detailed Implementation

[0029] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0030] Please see Figure 1 This invention provides a method for predicting the equivalence of steady-state pharmacokinetic extrapolation in complex long-acting formulations, comprising:

[0031] Step 1: Standardize the concentration-time, dosing regimen, and covariates of small samples with single or few doses, generate a standardized time, event identifier, and covariate dictionary, set censoring weights and population alignment weights, and form traceable modeling inputs and reference entries as the starting point.

[0032] The raw records of complex long-acting formulations often span different systems, sampling granularities, and units of measurement. Directly using these records for modeling can lead to time alignment errors, event ambiguities, and inconsistent covariate definitions, subsequently introducing systematic biases into later steps. Therefore, we first extract and unify time and events onto a single time axis, then map covariates to a set of computable and traceable standardized quantities. This process is repeated sequentially, completing time axis unification, event sequence encoding, and covariate dictionary standardization and mapping. To ensure comparability of all records on the same geometric timeline, we first construct a standardized timeline, then encode dosing and sampling events on this timeline, generating event identifier codes.

[0033] Furthermore, all subsequent kinetic equations can directly reference this unified coordinate system, avoiding relative temporal errors across files. Concentration records and dosing records have a mathematical one-to-one mapping, the time points for reaching stability and the recommended sampling window can be accurately replayed, and any covariate stratification can be reproduced on a unified time scale. Specifically, the standardized time is defined as follows:

[0034] Wherein: Subject number , No. Subjects; Sample point number The subject's first The sequence number of each time point; sampling time Subject number The absolute time of the first sampling; the time of the first dose administration. The absolute time of the subject's first administration; the length of the dosing interval. The subject's planned dosing interval; standardized time Dimensionless time coordinates, with a theoretical lower bound that can be less than zero (used to retrospectively trace the baseline before drug administration), and an upper bound determined by the study period.

[0035] Furthermore, in standardized time Generate event flag codes on the axis ,in Indicates a drug administration event. Indicates a sampling event. This indicates an invalid or revoked record; the encoding rules are written into the data processing record, and all revoked records retain their original traces and have the reason for revocation appended at the end.

[0036] In application, once the timeframe is unified, the concentration and dosing records of all subjects can be overlaid and viewed on the same coordinate system; after event coding, subsequent modeling can be performed based on event identifier codes. Directly identify the boundary between discrete dosing and continuous observation, reducing manual table matching; in the stability assessment stage, use standardized time... Using this as a benchmark, the relative position of the number of doses to reach stability can be directly located, improving the feasibility of the recommended sampling window.

[0037] Since covariates come from a wide range of sources and have different dimensions, they must be transcribed into a standardized covariate dictionary—a standardized covariate vector—with clearly defined units and value ranges. To avoid instability caused by common linear scaling mismatches in subsequent parameter estimation, a power-logarithmic mapping is used for key covariates, allowing variables of different orders of magnitude to enter a computable space of the same order of magnitude.

[0038] Therefore, covariates are mathematically comparable and additive. Covariate selection and stratification are carried out under a unified standard, and reference model items can be compared according to the same standard, reducing structural selection bias caused by different dimensions.

[0039] For each subject Raw values ​​such as weight, liver function score, kidney function score, and genetic typing are collected to form a covariate dictionary, recording its source, unit, and value boundaries. Discrete variables are stored in the database using deterministic encoding methods (e.g., genetic typing is represented by integer codes with an enumeration table). Power-logarithmic mapping is used for weight, and a standardized weight quantity is defined.

[0040] Among them: standardized body weight Dimensionless real number; numerical weight , No. The weight of the subjects; reference weight A pre-defined reference constant is used; liver and kidney function scores are scored using a pre-constructed scoring scale, and the scale version number is recorded simultaneously to avoid confusion; a standardized covariate vector is formed. .

[0041] In application, standardized covariates are brought to the same order of magnitude, reducing instability in subsequent numerical solutions; the one-to-one mapping between the covariate dictionary and the standardized covariate vectors ensures that any backtracking can trace back to the original record; reference model entries can be compared with... A point-by-point comparison was conducted under the same caliber.

[0042] Furthermore, early data for complex long-acting formulations often contain observations below the lower limit of quantitation. Simply deleting these observations would disrupt the continuity of the time series and exposure. At the same time, the distribution of sample covariates often deviates from the target population, and failure to correct this would lead to bias in steady-state exposure and equivalence.

[0043] Therefore, we first use left censoring weights to label the lower limit records and explain the numerical substitution strategy, and then use population alignment weights to correct the density ratio of the difference between the sample and the target population at the covariate level. By sequentially completing the labeling of the lower limit of quantitation and the generation of censoring weights, the construction of reference model entries and the estimation of population alignment weights, we can not destroy the factual records, but also provide a stable numerical entry point for subsequent estimation, and make the sample close to the target population in terms of covariates.

[0044] To reflect the measurement boundary without deleting records, lower limit observations are uniformly labeled and censoring weights are provided for subsequent estimations to control the influence of these points in the objective function.

[0045] Therefore, the lower limit information is preserved, the continuity is not interrupted, subsequent estimates are numerically more stable, and the sensitivity of steady-state extrapolation to the low concentration range is reduced. For each observation point, if the observed concentration is less than the lower limit of quantitation, a lower limit marker is assigned.

[0046] Censorship weights are defined as follows:

[0047] Wherein: Censorship weight The value can be either 0 or 1; the observed concentration , for the first The first subject Concentration observed in the first instance; lower limit of quantitation , No. The quantitative threshold reported by the laboratory in the batch of subjects, in units of [number] quantitation thresholds. Indicator operator Take when the conditions are met Otherwise, it is set to 0, which is used to apply special processing to the left censored points in the objective function of subsequent modeling.

[0048] In application, weights are removed. The explicit existence of the left censoring point allows for customized likelihood processing instead of deletion, conveying the lower bound fact while avoiding artificial breaks; the combination of lower bound marker and censoring weights makes it possible to interpret the source and processing of any playback.

[0049] Reference model entries are used to compare candidate structures, typical parameter ranges, and covariate effects under the same caliber. At the population level, to ensure that the distribution of sample covariates closely resembles that of the target population, population alignment weights are constructed using the density ratio approach. Therefore, reference entries provide anchors for structure selection, weight correction ensures that the samples represent the target population in a covariate sense, the bias in subsequent parameter estimation and extrapolation is suppressed, and the scenario simulation of steady-state indices is comparable at the covariate level.

[0050] Publicly available data were retrieved based on route of administration and active ingredient. Structural types, absorption / elimination characteristics, and covariate effect function forms were extracted and organized under the same criteria as the covariate dictionary—standardized covariate vector—to form a reference model entry list, ensuring that all covariates are shared. Each component defines the audience alignment weight:

[0051] Wherein: covariate components , No. The first subject The standardized values ​​of each covariate, corresponding to the standardized covariate vector. The components; the number of components It is a positive integer. Crowd alignment weight. , No. The dimensionless non-negative real number of the subjects;

[0052] Target distribution density , No. The marginal density function of each covariate in the target population can be extrapolated by the kernel density estimator. It is a non-negative real number and integrable over the entire domain.

[0053] Sample distribution density , No. The marginal density function of each covariate in the current sample; where: the marginal density is estimated using Gaussian kernel density; the bandwidth is minimized using the leave-one-out criterion; and crowd alignment weights are applied. Apply clipping:

[0054]

[0055] source: From step one; the target and sample observations are derived from the reference population and the current sample.

[0056] in: : The index of the covariate in the covariate dictionary, taken as... ; :No. A scalar value for each covariate (dimensionless, as standardization was completed in step one); The target audience is in the first... Marginal densities on each covariate are estimated using reference data from the target population. The current sample is at the . The marginal densities on each covariate are estimated using the sample data from this project.

[0057] When applying this approach, the reference entries align external knowledge with the project's standards to prevent structural assumptions from becoming detached from the data. Used to correct the density ratio of samples at the covariate level, making subsequent estimates and simulations closer to the expected target population; when combined with censoring weights When both covariates exist, they control both the measurement boundary and the representativeness of the covariates.

[0058] Step 2: Based on the modeling input, establish a population pharmacokinetic model that includes compartment structure, absorption and elimination mechanisms. Use weighted estimation and set a uniform quality threshold and covariate screening criteria to output an extrapolable and auditable final model and parameters.

[0059] Complex long-acting formulations exhibit absorption lag, release pathway superposition, and in vivo nonlinearity. If described by a single linear absorption, systematic bias is likely to occur in the early and late stages. At the same time, the distribution of sample covariates often deviates from the target population, and observations below the lower limit of quantitation are unavoidable. Direct deletion of points and unweighted estimation will both introduce extrapolation bias.

[0060] Therefore, the approach involves a two-stage linkage: first, identifying structures with dual-channel absorption and delayed expression capabilities based on a list of reference model items; and then aligning the population with weights. Deleted weights The likelihood is embedded to form an interpretable and extrapolable parameter estimation chain. This ultimately yields a final model structure and preliminary parameter range that are engineering-workable and statistically traceable, laying the foundation for subsequent covariate mechanism embedding and quality threshold validation.

[0061] To realistically depict the superimposed effect of donor release and receptor absorption in complex long-acting formulations, a dual-channel parallel absorption plus absorption lag description is adopted: one channel uses zero-order release to characterize drug release, and the other channel uses first-order absorption to characterize site absorption, with the actual onset of action delayed by the lag time.

[0062] Therefore, the structure can simultaneously accommodate sustained release and absorption. The early plateau and late tail segments are characterized simultaneously under a set of parameters, and the inter-cycle accumulation of dosing intervals can be reproduced within the same kinetic framework, providing a stable baseline for calculating the steady-state threshold and extracting the steady-state interval. Within the limits allowed by the reference model item list, it is determined whether the central and peripheral compartments coexist, and the distribution volume of the central compartment is fixed. Dimensional mapping and clearance rate The unit caliber.

[0063] Wherein, the absorption input rate is defined. :

[0064] Wherein: absorption input rate The instantaneous input rate entering the central chamber is a non-negative real number; the zero-level channel ratio A dimensionless parameter ranging from 0 to 1, representing the proportion of drug release channels; zero-order release rate constant. Positive real number, unit is quantity per unit time; indicator operator When time Not less than the absorption hysteresis time Set to 1 if the condition is met, otherwise set to 0; absorption hysteresis time Non-negative real numbers; proportion of primary channels , is a dimensionless parameter ranging from 0 to 1, satisfying ;

[0065] First-order absorption rate constant A positive real number, with the unit being the reciprocal of time; the amount of drug to be absorbed. The amount of drug at a location that changes over time, expressed in units of amount, satisfies the conservation of the body side.

[0066] In application, this driver unifies hysteresis, zero-order, and first-order into a single formula, which is convenient for numerical implementation and compatible with formulation mechanisms. It also provides intuitive interpretability when comparing candidate structures. Simultaneous fitting of long plateau and long-tail data segments is more robust, improving the stability of subsequent steady-state extrapolation curves.

[0067] Furthermore, to ensure that the sample representativeness and the measurement boundary are treated with the same caliber during the estimation stage, population alignment weights are applied. With censorship weight Simultaneously, the objective function is written, and for observations below the lower limit of quantitation, left censoring is used instead of direct deletion. Therefore, the distribution bias of covariates is suppressed, the information of the lower limit segment is preserved, the sensitivity of parameter estimation to the early low concentration segment decreases, and a parameter solution that is closer to the target population is formed.

[0068] For each subject Assigning weights to the sample layer For each observation point Give deletion instructions Construct weighted negative log-likelihood :

[0069] Wherein: weighted negative log-likelihood Non-negative real numbers; crowd alignment weights , No. Dimensionless nonnegative real numbers of subjects; censored weights A binary indicator, taking either 0 or 1; observed concentration. , No. The first subject Quantitative values ​​from each observation; individual parameter vectors , No. A set of parameters for the subjects, including clearance rate. Distribution volume First-order absorption rate constant Zero-order release rate constant Absorption lag time wait;

[0070] Left deletion of likelihood term When the observation is below the lower limit of quantitation Log-likelihood at time, For the first The lower limit of quantitation for each subject corresponds to the laboratory's limit of quantitation; where:

[0071]

[0072] Observed concentration : No. The number of subjects in the first Measured drug concentrations at specific time points; unit is concentration (e.g., ng / mL); range (If it is recorded as being below the lower limit of quantification, it will not be included in this formula, but will be deleted from the left side);

[0073] Model predicted concentration The final model determined in step two is applied to the subject. At the point of time Predicted concentration; unit same Source: Obtained by solving the differential equations of the central / peripheral chambers and superimposing the absorption inputs;

[0074] Individual parameter vector Subjects The parameter set must include at least the clearance rate. Central room distribution volume First-order absorption rate Zero-level release duration Absorption lag time In the case of a two-compartment structure, the intercompartmental clearance rate is also included. Peripheral chamber distribution volume Range: All the above kinetic quantities are positive; Unit: Use the same standard as specified in the instruction manual (e.g., L / h, L, ...). h);

[0075] Scale parameters The standard deviation of observation residuals on a logarithmic scale describes the intensity of observation noise; range: Dimensionless; Source: Sui It is estimated together in the weighted likelihood in step two;

[0076] Natural logarithm :by The base is used to convert multiplicative errors into additive errors, and the scope and units have been eliminated by logarithms;

[0077] constant : Mathematical constant, used as a regularization constant term for log-normal error. Index For the subject's serial number, The time point number; range: positive integer; one-to-one correspondence with data records. Process position explanation: This formula only sums for "unmarked and censored observation points" and proceeds to the weighted negative log-likelihood objective function in step two, used for estimation. and .

[0078] Left-censored log-likelihood, used for observations below the lower limit of quantitation:

[0079] Lower limit of quantitation Subjects The lower limit of quantitation for the experimental batch; units are the same. ;scope: Source: Laboratory report; if the batch level is consistent, it can also be recorded as a batch-level constant; model predicted concentration. Individual parameter vector Scale parameters The meaning is the same as above, taken from the same censored observation. Predictions and parameters. Standard normal distribution function. : Give the observation value under logarithmic scaling error not exceeding Cumulative probability; calling the standard normal cumulative distribution interface of commonly used mathematical libraries; indexing. Each observation point is assigned a unique value to the one marked as being below the quantitative lower limit; this value is then included in the summation using censoring weights in the objective function.

[0080] When an observation point is assigned a censoring weight in step one, the observation value is no longer calculated using the log-likelihood of the observation point, but instead uses the weighted negative log-likelihood objective function in step two.

[0081] Penalty coefficient Non-negative real numbers; regularity term For the full parameter set The smoothing regularization function, in its specific form, employs piecewise differentiable entropy constraints for stable estimation.

[0082] To stabilize the estimation and suppress extreme parameters, the regularization term employs entropy smoothing:

[0083] Wherein: the complete parameter set The set of all fixed effects parameters to be estimated; Each item is a scalar; weight Take positive real numbers (default) to This piecewise differentiability avoids sharp, undifferentiable points.

[0084] In application, this objective function combines representativeness correction, lower bound consistency, and regularization constraints into one, directly corresponding to engineering facts at the weight and likelihood level, avoiding sample pruning at the data level; it can converge stably in scenarios with a mixture of long tails and lower bounds, providing a reliable foundation for subsequent covariate mechanism embedding.

[0085] The sources of variation in complex long-acting formulations are extensive. If parameters such as weight, liver and kidney function, and genetic typing are not explicitly embedded, inter-individual differences will be misassigned to random terms, leading to systematic biases in terms of time to stability and steady-state exposure. At the same time, single graphical diagnosis is insufficient to cover biases at the temporal and quantile levels.

[0086] Therefore, the approach proceeds with covariate mechanism functionalization followed by triple threshold verification: first, standardized covariates are embedded in a log-linear form into the clearance rate and absorption-related parameters; then, robustness is verified using quantile checks and resampling, providing explicit, actionable pass / fail criteria. The final model achieves clear decoupling between parameters, covariates, and random terms, and meets the coverage requirements at the quantile level.

[0087] To allow the influence of covariates to enter the interpretable parameter channel, body weight was standardized using a log-linear function. The system embeds clearance rates and absorption terms into liver function scores, kidney function scores, and genetic typing codes, while preserving the exponential structure of inter-individual differences. Therefore, the influence of covariates is no longer mistakenly absorbed into the random term, the posterior distribution of individual parameters converges faster, and the steady-state extrapolation for different covariate stratifications provides a clear interpretation. The recommended number of times to reach stability and the sampling window are more transferable. For example, the individual clearance rate function is defined as follows:

[0088] Among them: individual clearance rate , No. Clearance rate per subject, in volume per unit time; typical clearance rate Typical values ​​for the population, in units of the same Standardized weight measurement Dimensionless real number, derived from the definition in step one; liver function score This is a standardized score, a standardized quantity of discrete or continuous scores, with dimensionless units and a range of values ​​defined according to a pre-defined scale; renal function score. For standardized scoring, a dimensionless standardized quantity; body mass index. Liver function coefficient Renal function coefficient The dimensionless coefficient to be estimated; inter-individual difference term. A log-scaled random term with zero mean is used to characterize unexplained individual differences.

[0089] Solving for the coefficients and variances of the random term using the same weighted objective function framework allows for coupled estimation of covariates, parameters, and random terms within a single solution round. In application, this function decouples the influence of covariates from the random term and places it onto interpretable coefficients, giving stratified steady-state predictions of different covariates a clear direction; when weighted with population data... When used together, representativeness and interpretability can be controlled simultaneously at the sample layer and the model layer.

[0090] Traditional methods that only consider mean deviation can easily mask long-tail or low-quantile deviations. Introducing a quantile coverage difference measure in addition to the prediction and verification plot allows for a weighted comparison of the full distribution under time conditions.

[0091] Therefore, the diagnosis is extended from point estimation to the distribution layer, the long plateau and long tail bias unique to long-acting formulations are explicitly exposed, the quality threshold no longer depends on a single indicator, and whether to enter the simulation has a clear pass / fail judgment.

[0092] Furthermore, the time-weighted quantile coverage difference is defined. :

[0093] Among them: time-weighted quantile coverage difference Non-negative real numbers; time weighting function It is a nonnegative integrable function used to emphasize critical time windows (such as plateau and tail segments), and its specific shape is given in the data processing record; conditional empirical distribution. At that moment Below, the concentration of the simulated sample is no higher than Empirical distribution function value; dimensionless; calculated using time kernel smoothing with uniform bandwidth (Gaussian kernel, bandwidth) ),and The same time grid and bandwidth are used.

[0094] Conditional empirical distribution At that moment Below, the concentration of the observed sample is not higher than Empirical distribution function value; dimensionless; calculation scope and Completely consistent; points marked as below the quantitative lower limit are not directly used for this distribution, but are processed according to the aforementioned censoring criteria and recorded in the data records.

[0095] Time weighting function : For any given moment within the study time window The assigned non-negative weights; used for the time-weighted integral in the quality threshold, dimensionless. Study the lower bound of the time window. Study the upper bound of the time window The start and end dates of the entire evaluation period; unit: time; meeting the following requirements. Dosing interval A fixed interval between two consecutive doses; the unit is time. The start time of the cycle. No. The starting point of each dosing cycle (usually the time when the drug is administered in that cycle); derived from an event-time series; unit is time; Number of cycles The number of complete or partial dosing cycles within the study time window; a positive integer. Plateau segment weight constant. : The weight assigned to the plateau segment; dimensionless and non-negative; generally used to emphasize the plateau region after stability is achieved. Tail segment weight constant. : Assign weights to the tail segment (end of the interval); dimensionless, non-negative; usually not less than Basic weight constant Weights for all time periods except the platform and tail segments; dimensionless and non-negative; typically satisfying the following conditions. Platform segment start and end coefficients Determine the relative position of the plateau segment in each dosing cycle to meet the following requirements: Dimensionless; plateau segment interval is , Tail segment starting point coefficient Determine the relative position of the start point of the tail segment in each dosing cycle to satisfy... and The tail segment interval is Union symbol For each dosing cycle within the study time window, the corresponding intervals are taken and then the union is performed to ensure complete coverage of the weight definition across multiple cycles. Difference set This refers to the remaining time interval after excluding the platform segment and the tail segment, which is uniformly assigned a value. .

[0096] Conditional distribution function In time Empirical distribution function of simulated concentration; conditional distribution function In time The empirical distribution function of the observed concentration; the integral domain covers the study time window and the measurable concentration domain.

[0097] Wherein: query time At any point within the study time window, measured in time; used to assess the concentration distribution in the vicinity of that moment. Observation time Sample No. Each actual sampling time is represented in time; each corresponds one-to-one with a single observation. Observed concentration. :exist The concentration measured at the point should be consistent with the data in terms of units; if it is marked as being below the lower limit of quantitation, that point will not be used in this formula and will be subject to censoring.

[0098] Threshold concentration The horizontal axis variable used in calculating the empirical distribution, taking values ​​ranging from the joint minimum to the maximum of the sample and simulated concentrations. Indicator function. : 1 is used when the observed concentration is not higher than the threshold, otherwise 0 is used. Time kernel Gaussian kernel, for distance Perform time-neighborhood weighting. Bandwidth The time smoothing scale is a positive real number. It can be minimized using the leave-one-out method, or set to a fixed value. The actual value is recorded in the data processing record. The result is calculated using measured data. The final model was calculated using simulation samples within the same time window and grid. .

[0099] Time-weighted quantile coverage differences Together with the goodness-of-fit plot for intuitive verification of bootstrap re-estimation stability, this serves as a triple threshold for entering the simulation. Failure to pass any one of these thresholds results in retrying at the structural candidate or covariate setting level. In application, this metric integrates time importance and distributional differences into a computable scalar, facilitating automatic determination within the engineering process. Combined with graphical diagnostics and bootstrap re-estimation, it can form robust quality thresholds at the quantile, graphical, and resampling levels.

[0100] Step 3: Call the final model and covariate effects, generate a virtual population according to the population covariate distribution, and conduct multiple dosing simulations; determine the stabilization by the threshold of relative change in trough concentration before two consecutive dosings, and determine the peak and trough sampling time window in the steady state interval based on the rate of change threshold.

[0101] The steady-state behavior of complex long-acting formulations is sensitive to the distribution of covariates, dosing intervals, and absorption hysteresis. If the virtual population fails to maintain the correlation structure among covariates, or if dosing events do not strictly fall on a unified coordinate system, subsequent determinations of stability will exhibit systematic biases. Therefore, we first ensure consistency with the target population and preserve the relevant structure at the covariate level, and then align dosing events with the computational grid at the time level to ensure that the absorption input and solution step size are consistent. This ultimately yields a multiple-dose curve with clearly defined covariate definitions and a unified time axis, forming a unique baseline for subsequent determinations.

[0102] To avoid structural breaks caused by sampling solely based on marginal distributions, the covariate vector is reconstructed in the standardized space according to the covariate correlation structure of the target population, and then mapped to the typical value channels of individual parameters.

[0103] First, the reconstruction of the relevant structure is given, so the linear correlation between covariates can be restored in a mathematical sense, thereby ensuring that the consistent input to the parameter channel at the individual level is generated by the covariate vector entering the parameter channel through the function established in the previous step to generate a virtual population parameter set that can be called by the solver.

[0104] Let the mean vector of the dummy covariates be... The lower triangular decomposition matrix is The standard normal vector is ,structure:

[0105] Where: dummy covariate vector , No. A standardized covariate column vector for each virtual individual, with element dimensions consistent with the covariate dictionary from step one; a virtual covariate mean vector. The target mean column vector of the covariate, with each element being a real number; lower triangular decomposition matrix. The Joliski decomposition of the target correlation matrix, a lower triangular matrix with real elements; standard normal vector. An independent standard normal column vector of the same order as the covariate dimension, with each element taking values ​​that follow a distribution with a mean of 0 and a scale of 1, used to introduce relevant structures.

[0106] Furthermore, the virtual covariate vector Substitute each item into the parameter-covariate function (e.g., individual clearance rate function and absorption-related function) specified in step two, and combine with the parameter vector. Covariate coefficient set Set of variances among individuals The parameter set for each virtual individual is generated according to the exponential structure and stored on disk as a list of virtual individual parameters.

[0107] In application, virtual covariate vectors are used. The construction of the covariate correlation is maintained in the virtual population, avoiding parameter mapping distortion caused by marginal sampling alone; the parameter-mapped population has the same input caliber as the target population, providing a stable entry point for curve generation.

[0108] To ensure that the absorption drive is consistent with the solution step size, the dosing events are generated with precise timestamps and these timestamps are aligned with normalized time-event flag codes. The event sequence is then fed into the ordinary differential equation solver along with the absorption input function.

[0109] Among them, the first dose is generated based on the absolute time of the first dose. The absolute time of the first dose:

[0110] Wherein: Dosing time , No. Within the first dosing cycle The absolute moment of each actual dosing event; the start time. The absolute time of the first dose; the dosing interval. Fixed positive real number, unit is time; offset The non-negative time offset of multiple events within the same period ensures that multiple events are within the same interval and do not overlap.

[0111] Furthermore, the timing of drug administration... With standardized time Event flag code Align, generate event-time series, and input the absorption function from step two. By binding this sequence, a variable-step double-precision solver for ordinary differential equations is invoked to integrate to the predetermined end time of the study, obtaining continuous concentration-time curves for multiple dosing administrations. In application, the explicit expression of event times ensures strict consistency between the dosing rhythm and the solution grid, and absorption lag and multi-channel release can be accurately expressed numerically; , Alignment ensures consistency of caliber across steps.

[0112] The rate at which complex, long-acting formulations reach steady-state conditions is influenced by clearance, absorption lag, and dosing rhythm, making empirical methods using fixed dosing frequencies prone to misjudgment. Furthermore, the slope of the curve near the peak and trough values ​​determines the sampling tolerance; without setting windows based on the rate of change, peak or trough values ​​may shift. Therefore, we first identify the trough value in each dosing cycle and determine the number of times steady-state conditions are reached using a relative rate of change threshold. Then, within the steady-state range, we define sampling time windows for peak and trough values ​​using the rate of change threshold. This ultimately results in a list of peak and trough sampling windows for the distribution of dosing frequencies at steady-state conditions, providing a definite time window for steady-state non-compartmental analysis and equivalence determination.

[0113] At the end of each dosing interval, the trough concentration sequence before the next dosing is extracted, and the relative change rate between two adjacent troughs is used as the criterion. Stabilization is considered achieved when the relative change rate is below a threshold. The specific criterion is as follows:

[0114] Among them: the first Valley concentration before the first dose In the The concentration at the end of the dosing interval and before the start of the next dosing; Valley concentration before the first dose Definitions and units are the same as before; threshold Dimensionless positive real number, default value is taken as It allows for limited adjustments to the formulation and design within a pre-defined allowable range.

[0115] Furthermore, on the multiple dosing curves for each virtual individual, according to event time... The concentration at the end of each period is extracted as the trough sequence; and, starting from the second period, the relative rate of change is calculated for each term and compared with a threshold. For comparison, the first time the condition is met is recorded as the number of times the stable condition is reached. In application, the relative rate of change is used instead of the absolute difference to avoid interference from baseline differences between individuals in the determination; the criterion of stopping the drug upon first meeting the condition is easy to implement in engineering and can give a clear number of times the stable condition is reached for each virtual individual.

[0116] To provide an executable peak and trough sampling time window within the steady-state range, the time window boundary is defined based on the local rate of change of the curve: when the absolute value of the rate of change is not higher than a fixed proportion of the reference slope, the representativeness of that time period for the peak or trough is considered acceptable. Wherein, the following definition applies:

[0117] Wherein: sampling time window set The set of times that satisfy the conditions; lower bound of time. With the upper bound of time The start and end times of the steady-state interval; steady-state concentration function. The concentration-time function over a complete dosing interval after reaching steady state, expressed in units of concentration; specifically, the concentration-time function over any dosing interval after reaching steady state. The steady-state concentration function is: And satisfy the periodic extension condition:

[0118] The state vector is: The quantities are, in order, the dosage at the location, the dosage in the central chamber, and the dosage in the peripheral chamber; This represents the central chamber drug quantity component in the steady-state one-period solution. System matrix and periodic input:

[0119] ;

[0120] Analytical expression for the initial value of the periodic boundary:

[0121] ,in Let the identity matrix be '[');'. Substituting this into the above equation yields the unique periodic steady-state solution, which is then... Obtain the concentration trajectory.

[0122] Let be the steady-state vector, and let the matrix exponent be . Implemented using common numerical linear algebra methods; convolution integrals are implemented using Gaussian or trapezoidal quadrature. In discrete implementation, it is equivalent to... to Instantaneous increase ; For constant injection within the interval, it is treated as a constant source term during numerical integration. Phase is adjusted during periodic extrapolation. ,Pick .

[0123] steady-state concentration function Central chamber concentration at any point after stabilization; unit consistent with experimental concentration; domain is real time; local dosage. Central room dosage Peripheral chamber dosage : Quantity of medicine in each room; unit is quantity; represents the state vector components; system matrix. : A constant coefficient matrix consisting of the clearing and distribution dynamic constants; dimension is ; Periodic input : Input vector within a dosing interval; the first component is the pulsed input of the first-order channel, the second component is the constant input of the zero-order channel, and the third component is zero; first-order absorption rate : Reciprocal of time; estimated by the model; clearance rate Volume per unit time; from model estimation; intercompartmental clearance rate Volume per unit time; represents the exchange intensity between the central chamber and the peripheral chambers; central chamber distribution volume. Peripheral chamber distribution volume Volume; Zero-order scale First-level proportion Dimensionless; the sum of the ratios of the two channels is 1; dosage : Amount administered per dose; unit is dose; given by the dosing regimen or set by simulation; absorption lag time Time; Inflow begins after lag; Duration of level zero Time; length of the zero-order constant intake interval; dosing interval Time; fixed length of adjacent dosing intervals; starting point of the steady-state interval. Time; the starting point of the interval determined by the stability assessment in step three; the initial periodic state. The steady-state one-period solution is in The state vector; analytically obtained from the above equation; the identity matrix. :and Identity matrices of the same order; Dirac pulse In analytical expressions, this represents instantaneous injection; in numerical implementation, it represents instantaneous increment operations; it is an indicator function. : A function that takes 1 within an interval and 0 outside the interval; used to represent the time window of zero-level input.

[0124] proportionality coefficient A dimensionless positive real number is used to adjust the sampling window width, with a value between 0 and 1; the reference slope ,exist The slope represents a non-negative real number, specifically within that interval. The median is used to stabilize the sensitivity to abnormal spikes.

[0125] Specifically, numerical differentiation is performed on a complete interval after the number of doses administered to obtain the rate of change sequence, and then the calculation is performed. Furthermore, the sampling windows are expanded to both sides around the peak and valley times, respectively, and continuous intervals that satisfy the set definition are taken as peak sampling windows and valley sampling windows, and a list of sampling windows is output.

[0126] In application, the rate of change threshold transforms representativeness into a computable set criterion, avoiding the execution difficulties caused by excessively narrow windows near peaks and valleys, as well as the representativeness reduction caused by excessively wide windows; the median reference slope suppresses the influence of abnormal peaks, and the sampling window is more practical in engineering.

[0127] Step 4: Perform non-compartmental analysis within the complete dosing interval after stabilization to calculate the dosing interval area, steady-state peak value, and steady-state trough value; conduct interval determination and efficacy-sample size linkage simulation for the steady-state index of the test and reference, and output a consistent research design and report.

[0128] The curve morphology (plateau-peak-trough) of complex long-acting formulations within one dosing interval after reaching stability determines the exposure level and the reliability of representative sampling. If the steady-state interval and peak-trough sampling window are not fixed at the document level, the transitional information before stabilization will be mistakenly included in the statistics, thus affecting the equivalence judgment.

[0129] Therefore, the steady-state threshold and the upper limit of the dosing interval are first locked on the time axis, and peak-valley sampling is constrained within this interval using a sampling window list. Then, the dosing interval area is calculated using a trapezoidal integral method (non-compartmental analysis), and steady-state indices corresponding to the peak-valley positions are extracted simultaneously. This results in a steady-state indices list, with all entries including the number of times steady-state was reached and the time window source trajectory, allowing direct entry into equivalence statistics.

[0130] To ensure consistency, the starting point of the steady-state interval is determined by the number of dosing cycles to reach stability in step three, and the first complete dosing interval after this starting point is recorded as the lower time bound. With the upper bound of time Therefore, any subsequent indicators are in Internal calculations are performed to avoid interference from the transition period on exposure estimates. Peak-valley sampling is only performed within the time window allowed by the sampling window list, and all steady-state indicators have replayable time sources.

[0131] This includes reading the stability assessment report, locating the number of doses administered to each virtual individual to determine stability. and ; and, in The sampling grid is generated with the sampling window list as the priority constraint. When additional computational nodes exist outside the window, it is only used for integral interpolation and not for peak-valley statistics. In application, the time boundary and sampling grid are fixed in the same file to ensure that no caliber drift occurs when the next stage is read in; the peak-valley positions are constrained by the sampling window list to reduce execution risks.

[0132] In the locked The area within the dosing interval is calculated using the trapezoidal integral method, and steady-state indices corresponding to peak and trough values ​​are extracted according to the sampling window list. The dosing interval area is defined as follows:

[0133] Wherein: dosing interval area , for in The area within, expressed as concentration multiplied by time; steady-state concentration function Within the steady-state interval, at time The concentration value, in units of concentration; sampling time Time nodes in the sampling grid, in units of time; number of grid points Positive integer; dosing interval Fixed positive real numbers;

[0134] Within the peak and trough time windows defined in the sampling window list, the maximum concentration is searched as the steady-state peak concentration, and the minimum concentration is searched as the steady-state trough concentration, with corresponding timestamps recorded for playback. In application, area calculation and peak-trough extraction share the same time metric to avoid overlapping metric; both peaks and troughs are bound to timestamps for easy item-level backtracking.

[0135] The steady-state equivalence of complex long-acting formulations depends on both the ratio of exposure indicators and the design's ability to accommodate inter-individual variability. Without linking equivalence statistics with efficacy-sample size, feasible design conclusions cannot be obtained before project initiation. Therefore, this approach first calculates the logarithmic mean difference and geometric mean ratio for steady-state indicators at the item level, forming confidence intervals. Then, repeated sampling is performed on a two-dimensional grid of inter-individual variability and sample size. The frequency of statistical intervals falling within equivalence intervals is used as efficacy values, and designs are recommended based on this. Finally, a steady-state equivalence report and efficacy-sample size curve are generated, along with recommended schemes with sampling windows and timestamp sources.

[0136] For each steady-state index (dose interval area, steady-state peak concentration, steady-state trough concentration), the logarithmic mean difference between the test formulation and the reference formulation was calculated, and the geometric mean ratio was determined. The equivalence conclusion was based on the confidence interval falling within the equivalence interval. Here, is defined as: Where: logarithmic mean difference The difference between the logarithmic means of the two formulations; the ratio of their geometric means. Equivalence evaluation ratio; mean of logarithmic steady-state index The test formulation in steady-state parameters Logarithmic mean of the above, Take the dosing interval area, steady-state peak concentration, and steady-state trough concentration; logarithmic steady-state index mean. The logarithmic mean of the reference formulation on the same indicator;

[0137] Furthermore, to estimate the standard error With degrees of freedom percentiles Calculate the 90% confidence interval:

[0138] Among them: 90% confidence interval Confidence interval for the geometric mean ratio; standard error of estimation Standard error of the logarithmic mean difference, in dimensionless units; percentiles The degrees of freedom are of Distribution in cumulative probability Quantiles at; degrees of freedom A positive integer, determined by the design and sample size. Judgment criterion: When the 90% confidence interval... All contained in the interval When the time comes, record the equivalent of this indicator.

[0139] In application, the logarithmic scale stabilizes the symmetry of the ratio, and the criteria for falling within the interval are consistent with regulatory practices; the item-level interval reporting and timestamp mapping enable subsequent reviews to be replayed at the specific individual-specific sampling window level.

[0140] To provide clear recommendations on sample size and design type during the project initiation phase, repeated sampling was performed on a two-dimensional enumeration grid of inter-individual differences and sample size, with the frequency of interval falls used as the efficacy value, and cost and ethical burden used as auxiliary screening conditions.

[0141] Wherein: Define efficacy value Monte Carlo estimates:

[0142] Among them: efficacy value The probability of equivalence success under a given design; the number of repetitions. A positive integer, representing the number of rounds of resampling; indicator function The value is 1 if the condition within the parentheses is true, otherwise it is 0; 90% confidence interval of rounds The interval calculated after generating samples in this round; the lower and upper bounds of the equivalent interval. Fixed as and .

[0143] Specifically, the calculation is performed on the combination of the enumerated inter-individual difference level and sample size. We screened combinations that met the efficacy threshold of at least 80%, and provided recommendations for crossover or parallel designs based on engineering constraints such as the total number of blood collections and the study period.

[0144] When applied, the frequency interpretation of efficacy values ​​is intuitive and clear, making it easy for management to understand; when combined with cost and ethical constraints, the output design parameters are naturally executable.

[0145] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0146] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0147] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or units may be electrical, mechanical, or other forms.

[0148] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0149] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for predicting the equivalence of steady-state pharmacokinetic extrapolation in complex long-acting formulations, characterized in that: include, Collect human concentration-time, dosing regimen and covariates for single / few doses and standardize the data. Generate censoring weights for low-limit observation labels and calculate population alignment weights. Output the modeling input dataset, data processing records and reference model entries. Based on the modeling input dataset, candidate structures including room number, absorption and elimination forms are set. Nonlinear mixed effects are used to combine two types of weighted estimation parameters, and the final model and covariate effects are determined according to the quality threshold. Based on the final model and covariate effects, a virtual population was generated according to the population covariate distribution and multiple dosing simulations were performed; the threshold of relative change in trough concentration before two adjacent dosings was used to determine the steady state, and the peak and trough sampling time windows were determined according to the rate of change threshold in the steady state interval; Within the complete dosing interval after reaching stability, the dosing interval area, steady-state peak value, and steady-state trough value are calculated. The geometric mean ratio of the steady-state indices of the test and reference is calculated and determined according to the confidence and equivalence intervals. The efficacy is simulated under the sample size and enumeration of inter-individual differences, and the recommended design is output. Among them, for each subject Assigning weights to the sample layer For each observation point Give deletion instructions Construct weighted negative log-likelihood : Among them: weighted negative log-likelihood Non-negative real numbers; crowd alignment weights , No. Dimensionless nonnegative real numbers of subjects; censored weights A binary indicator, taking either 0 or 1; observed concentration. , No. The first subject Quantitative values ​​from each observation; individual parameter vectors , No. A set of parameters for the subjects, including clearance rate. Distribution volume First-order absorption rate constant Zero-order release rate constant Absorption lag time ; Left deletion of likelihood term When the observation is below the lower limit of quantitation Log-likelihood at time, For the first The laboratory's lower limit of quantitation for each subject; penalty coefficient. Non-negative real numbers; regularity term For the full parameter set The smoothing regularization function, in its specific form, employs piecewise differentiable entropy constraints for stable estimation; To stabilize the estimation and suppress extreme parameters, the regularization term employs entropy smoothing: Where: the set of all parameters The set of all fixed effects parameters to be estimated; Each item is a scalar; weight Take positive real numbers; Segmented differentiability.

2. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 1, characterized in that: Establish a covariate dictionary and unify units and codes. Determine the population alignment weights based on the ratio of covariate marginal density between the target population and the sample, and perform interval pruning. Mark censoring weights for records below the lower limit of quantitation. Generate modeling input datasets and data processing records with version fingerprints.

3. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 2, characterized in that: Candidate structures include parallel zero-order release and first-order absorption with hysteresis time, and inter-chamber exchange between the central chamber and peripheral chambers; parameter estimation is performed using weighted likelihood with incorporation population alignment weights and censoring weights; observation errors are modeled on a logarithmic scale; and the final model and covariate effect functions are output.

4. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 3, characterized in that: The quality threshold includes at least the following: generating a goodness-of-fit plot, a time-stratified prediction validation plot, and bootstrap parameter reestimation; covariate screening is performed using forward and backward thresholds and a fixed random seed; and upper limits are set and recorded for the uncertainty of the main parameters to form an evidence package for passing the quality threshold.

5. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 1, characterized in that: Reconstruct the covariate correlation structure based on the covariate dictionary and generate a list of virtual population parameters; generate event-time series according to dosing intervals and start times and align them with standardized time and event flag codes; The random seed and event list are fixed by document logging.

6. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 5, characterized in that: The determination of stabilization is based on the relative change of trough concentration before two consecutive administrations being lower than a preset threshold. Within the complete dosing interval after stabilization, the sampling time window for peak and trough values ​​is determined based on the proportion of local change rate of the concentration curve not exceeding the reference slope, and a list of sampling windows with timestamps is generated.

7. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 1, characterized in that: Within the first complete dosing interval after stabilization, a sampling grid was constructed with the sampling window list as the priority constraint. The area within the dosing interval was calculated using trapezoidal integrals, and steady-state peak values ​​and steady-state valley values ​​were extracted. Peak values ​​and valley values ​​were searched and determined one by one within the corresponding time window, and the corresponding time information was recorded.

8. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 7, characterized in that: The geometric mean ratio and confidence interval of the steady-state indices of the test formulation and the reference formulation are calculated on a logarithmic scale, and the equivalence is determined by the inclusion of the preset equivalence interval; wherein, the degrees of freedom and design type are matched and fixed in the equivalence determination report.

9. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 8, characterized in that: Efficacy simulations were performed on a two-dimensional enumeration grid of inter-individual differences and sample size, with the number of repetitions and random seeds fixed. Among combinations that meet the efficacy threshold, engineering constraints of total blood collection and study period were superimposed to screen the study design, and the recommended design was output in document form.

10. The method for predicting the steady-state pharmacokinetic extrapolation equivalence of complex long-acting formulations as described in claim 9, characterized in that: The entire process uses the modeling input dataset, reference model item list, data processing records, stability determination report, sampling window list, steady-state index list, and equivalence determination report as the sole input and output criteria, and maintains a unique mapping between terms and variable identifiers, all of which are written into version fingerprints and timestamps.

Images