Methods and devices for personalized administration of pharmaceutical agents that enhance safe delivery within the therapeutic range - Patents.com

Abstract

A technique for generating an individualized dose protocol includes receiving first, second, and third data. The first data represents a continuous multivariate model having at least one distribution parameter characterizing population variation for dose-response to a drug. The second data represents a therapeutic value range for the dose-response. The third data represents a cost function based on the distance of the dose-response from the therapeutic range. The method includes evaluating expected costs for candidate doses based at least in part on a Koopman transformation of the distribution parameters and the cost function. When the expected costs are below a threshold cost, the candidate dose of the drug is administered to the subject.

Description

[Technical field]

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims priority to Provisional Patent Application No. 62 / 891,602, filed August 26, 2019, the entire contents of which are incorporated by reference herein under 35 U.S.C. § 119(e) as if fully set forth herein. [Background technology]

[0002] Different patients process administered drugs differently, and for some drugs, or some neonates, or some combinations, variability in delivery of drug concentrations at the same dose may be highly relevant to the efficacy of the drug. The term pharmaceutical drug as used herein means any substance administered to a subject for therapeutic effect, including drugs and other biological agents such as large molecules, nucleic acids, viruses, and bacteria. In general, it can be said that for some pharmaceutical drugs, careful and precise management of administration and delivery is very important. The process of managing the administration and delivery of drugs is called therapeutic drug management (TDM) in the healthcare industry. Here, the term TDM is used to refer to the management of both drugs and other pharmaceuticals. Summary of the Invention [Problem to be solved by the invention]

[0003] If a drug is delivered at too low a dose, the therapeutic effect may be compromised, whereas if it is administered at too high a dose, adverse and even dangerous side effects may occur. [Means for solving the problem]

[0004] Techniques have been provided for more efficiently determining the proper dosage of a drug during individualized treatment of a subject, which dosage increases the chances of avoiding costly deviations from the target therapeutic range.

[0005] In a first set of embodiments, a method automatically executed by a processor for generating an individualized dose protocol includes receiving first data indicative of a continuous multivariate model having at least one distribution parameter characterizing a population variation for a dose response to a pharmaceutical product. The method also includes receiving second data indicative of a therapeutic value range for the dose response. The method further includes receiving third data indicative of a cost function based on the distance of the dose response from the therapeutic range. In addition, the method includes evaluating expected costs for the candidate doses based at least in part on a Koopman transformation of the distribution parameter and the cost function. In addition, the method includes administering the candidate dose of the pharmaceutical product to the subject when the expected cost is below a threshold cost.

[0006] In various embodiments of the first set, evaluating the expected costs further includes performing a multi-dimensional quadrature integration.

[0007] In some embodiments of the first set, the at least one distribution parameter is based on a distribution of values ​​in an electronic health record (EHR).

[0008] In some embodiments of the first set, the drug is of a type selected from the group including antiepileptic drugs, anticoagulants, antiarrhythmic drugs, and immunosuppressants.

[0009] In some embodiments of the first set, the distance of the dose-response from the therapeutic range of the cost function is indicated by the therapeutic index.

[0010] In another set of embodiments, a computer-readable medium or system is configured to perform one or more steps of one or more of the methods described above.

[0011] Further aspects, features, and advantages will be readily apparent from the following detailed description, by way of a brief illustration of several specific embodiments and implementations, including the best mode contemplated for carrying out the invention. Other embodiments are capable of other and different features and advantages, and its several details can be modified in various obvious respects, all without departing from the spirit and scope of the present invention. Accordingly, the drawings and description should be regarded as illustrative in nature, and not as restrictive.

[0012] Embodiments are illustrated by way of example, and not by way of limitation, in the figures of the accompanying drawings in which like reference symbols refer to similar elements and in which: [Brief description of the drawings]

[0013] [Figure 1] FIG. 1 is a block diagram illustrating an exemplary pharmacokinetic and pharmacodynamic system model used in accordance with one embodiment. [Figure 2A] 1A-1C are a pair of graphs illustrating an example relationship between an initial condition probability density function (pdf) and a resulting dynamically variable probability density function (pdf), respectively, according to one embodiment. [Figure 2B] 1 is a set of plots illustrating the interaction of uncertainty in initial conditions, represented by pdf, with a cost function using a Koopman adjoint process, according to one embodiment. [Diagram 3] 1 is a flow diagram illustrating an exemplary method for objectively and automatically determining a dosing schedule, according to one embodiment. [Figure 4] 1 is a graph illustrating an example convergence of a Koopman Expectation calculation compared to a previous method using Monte Carlo Expectation, according to an example embodiment; [Diagram 5] FIG. 1 is a block diagram illustrating an example computer system upon which one embodiment of the present invention may be implemented. [Figure 6] FIG. 1 illustrates a chipset in which an embodiment of the present invention can be implemented. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0014] A method and apparatus are described for determining an individualized protocol for administering a pharmaceutical agent likely to remain within a therapeutic range in a subject during the treatment of the subject. In the following description, for purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the invention. However, it will be apparent to those skilled in the art that the invention may be practiced without these specific details. In other instances, well-known structures and devices are shown in block diagram form in order to avoid unnecessarily obscuring the invention.

[0015] Although the numerical ranges and parameters indicating broad ranges are approximations, the numerical values ​​set forth in the non-specific limiting examples are reported as precisely as possible. However, numerical values ​​inherently contain certain errors necessarily resulting from the standard deviation found in the respective test measurements at the time of writing. Moreover, unless otherwise clear from the context, numerical values ​​presented herein have an implied precision given by the least significant digit. Thus, a value of 1.1 means a value from 1.05 to 1.15. The term "about" is used to indicate a broader range centered on a given value, and unless otherwise clear from the context, it means a broader range around the least significant digit, for example, "about 1.1" means a range of 1.0 to 1.2. When the least significant digit is unclear, the term "about" refers to a factor of 2, for example, "about X" means a value in the range of 0.5X to 2X, for example, about 100 means a value in the range of 50 to 200. Moreover, all ranges disclosed herein are understood to encompass any and all subranges contained therein. For example, a range of "less than 10" for a positive-only parameter can include any and all subranges between (and including) a minimum value of 0 and a maximum value of 10, i.e., any and all subranges having a minimum value greater than or equal to zero and a maximum value less than or equal to 10, e.g., 1 to 4.

[0016] In various embodiments, the corresponding precision therapeutic toolset will be directly applicable to all medicines that benefit from TDM. For example, predictive models other than precise dosing can be incorporated to assess the appropriateness of palliative chemotherapy in metastatic cancer patients with poor prognosis by predicting survival, to use appropriate antidiabetic therapy in type 2 diabetes patients with different pathophysiology, or to predict the risk of cardiovascular events in a given period of time.

[0017] 1. Overview Biologically and pharmacologically plausible system models, such as pharmacokinetic (PK) / pharmacodynamic (PD) (PK / PD) and physiologically based pharmacokinetic (PBPK) models, are tailored to the individual subject dose concentration levels. System models for precision treatment are highly nonlinear and require multiple levels of uncertainty and variability, leading to the use of nonlinear mixed effects (NLME) models. NLME for statistical inference is demonstrated in a large literature of applied statistics. Furthermore, the form of the NLME model limits the type of patient data that can be used to easily provided quantities such as age, sex, and weight. To determine values ​​for the various categories of model parameters, the system models listed here may be used in conjunction with NLME in conjunction with modern machine learning algorithms such as deep (convolutional code) neural networks, density-based clustering of diffusion mapping data, or gradient boosting regression trees. This pushes us well beyond modeling of mean values ​​into a space filled with complex multidimensional healthcare information by both fixed and random effects models. Such modeling seeks to capture a therapeutically important subset of the underlying complex pharmacological mechanisms that can be used to individualize treatment one patient at a time. A subject as used herein is a human or animal patient, or a healthy subject or surrogate or simulation thereof.

[0018] 1.1 Background on nonlinear mixed-effects (NLME) models Figure 1 is a block diagram illustrating an exemplary pharmacokinetic and pharmacodynamic system model 100 used according to one embodiment. In Figure 1, rectangular nodes represent variable parameters, circles represent latent (unfilled) or observed (filled) random quantities, diamonds are deterministic given inputs, and unbounded nodes are constants. The system model 100 represents one or more observable drug levels y for individual i at observation time element j. ij Such a system is defined by three components: random effects, dynamics, and observations.

[0019] Symbol η i The individual random effect, denoted by η , captures the inter-individual variability (IIV) of individual i within the population of subjects. The variable η i represents one or more random effects for individual i, and is therefore typically a vector. The distribution of this effect across the population is parameterized by Ω. For example, as an individual patient visits the clinic on multiple occasions for tissue samples to obtain drug levels (such as blood draws to obtain drug concentration values ​​at that time point in the blood), К is the random effect of individual i on each occasion k. ik Other individual random effects may also be incorporated, such as the inter-opportunity variation (IOV) effect, represented by ik represents one or more between-chance random effects for individual i and is therefore generally a vector. The distribution of between-chance effects across the population is denoted by Π. The most common choice of random effect components is the formula η i =N(0,Ω) and K ik =N(0,Π) is some kind of multivariate normal distribution N with zero mean.

[0020] Distribution parameters as used herein characterize population variability regardless of origin, e.g., inter-individual variability, such as variability based on covariates, individual event-to-event variability, distribution around population means for pooled data, and measurement uncertainty. For example, after performing likelihood estimation in something like a Bayesian sense, η in Ω terms iWe learn the distribution of θ, but can also obtain the distribution of the θ term. Maximum likelihood estimation returns (approximately) the maximum value of this distribution, but Bayesian fitting methods do not just return a single point, but a point with uncertainty (similar to estimation with error bars). Therefore, we can introduce uncertainty into the θ term by incorporating those error bars. These errors can come from many sources, such as measurement errors in the data, numerical errors in nonlinear routines, etc. Also, in some embodiments, we do not make individual-specific estimates, but we want to do optimal dosing under uncertainty about the estimation errors in fitting θ.

[0021] Dynamics controls the time course of the system. ij is a series of time points t for individual i and time element j ij represents the unobserved (latent) dynamics of the model that evolves in u ij represents one or more latent dynamic variables of individual i and is therefore generally a vector. In some embodiments, these dynamics are represented by a vector Z i one or more non-random parameters of η i , and vector Z i The covariate variables Z are provided as the solution to a system of differential equations involving subject-specific zero or greater values ​​for the covariate variables. i is the modeled y ij is a vector of zero or more observable characteristics of individual i that are known to vary by η. In other embodiments, the dynamics are modeled in other ways, such as using continuous and discrete Markov chains, or stochastic, delayed, hybrid, and differential algebraic variants of differential equations. In some embodiments, the first dose estimate is based on η i There are no random effects represented by

[0022] The functional relationship between these variables is given by Equation 1a.

number

number

number

[0023] Observed value y ij is a vector of one or more actual observed response variables (also called dependent variables). These are determined by the observation function Ψ, which is ij Associated with u ij and the random error represented by the vector σ, whose distribution over the population is denoted by Σ, as given in Equation 1b.

number

[0024] The output of a nonlinear mixed effects (NLME) model is determined as follows: A random variable (η i ,K ik ,σ) are defined. Similarly to the parameter value θ for the population, the covariates Z i These parts then define the parameters of the dynamic model for the particular set of differential equations to be solved (analytically or numerically) as expressed in Eq. 1b. Then, the observation function Ψ is applied as expressed in Eq. 1c to obtain the error value

number

number

[0025] For training and individualized administration of medication, the symbol d ij =d i (t ij ) we are interested in the inverse problem of learning more about a model given real observed data. The goal is to i , and Ω, Π, Σ (hence effectively the amount of randomness in the random effects, random chances, and random measurements), or more explicitly, the distributions η i , K ik , and

number

number

[0026] In the PK compartment model, the drug concentration in a compartment (C) is equal to the amount (A) divided by the volume (V). For two compartments, C1=A1 / V1 and C2=A2 / V2. The drug concentration in the central compartment (C1) is equal to the concentration in plasma (CP). The clearance (CL, in liters L / hour h) is expressed by the fractional rate constant k a (units are per time, not to be confused with occurrence k). In pharmacokinetics, the volume of distribution is given in volume units (L) and the rate constant is the ratio of clearance to volume of distribution, e.g., k a For oral administration of a drug, at time t=0, the amount of drug in the central and peripheral compartments is zero (A1(0)=A2(0)=0) and the initial amount in the gastrointestinal tract (effective dose) is AGI(0)=D×S×F, where D is the administered dose of the drug, S is the salt factor (fraction of the administered dose that consists of pure drug), and F is the bioavailability coefficient (fraction of the dose that reaches the systemic circulation). Commonly used parameters in PK models include D=dose, τ=administration interval, CL=clearance, Vd=volume of distribution, k e = disappearance rate constant, k a= absorption rate constant, F = absorption rate (bioavailability), K0 = infusion rate, T = infusion time, and C = plasma concentration.

[0027] In an exemplary embodiment, one or more of these parameters may be explained or influenced by the distribution of typical values ​​of the random value for different individuals.

[0028] Electronic health records (EHRs) provide a rich data set for employing machine learning methods. Thus, in various embodiments, the numerical solution of the differential equations of the NLME model is implemented as the activation layer of a neural net architecture. By directly encoding the ODE dynamics and events into the neural net layer via its activation function, the properties of the modeling-based approach can be preserved so that the output values ​​(predicted concentrations) are biologically realistic and follow the proposed clinical model. Population parameters and random distribution functions are learned during this training while maintaining the ODE relationships. Then, in some embodiments, this dynamic layer is incorporated into an end-to-end neural network with individual treatment effect measurements as inputs and expected dose-response and its confidence (probability) as the output layer. In other embodiments, other methods of determining the model parameters θ are used.

[0029] The parameters of this system become the input nodes of the dynamic layer, which are chained from previous layers of the neural network as predictions from any data set. In this way, relationships between electronic health record (EHR) records, imaging results, and molecular biology assays (protein assays, SNP arrays, RNA sequencing) can be programmatically trained and estimated without the need for modeler assumptions, generalizing PK / PD modeling methods directly to modern big data.

[0030] For example, in a Bayesian estimation process, we consider the dynamical system u underlying multiple individuals i. i Observed value y iGiven the likelihood of a given fixed effect θ and the random effect variance Ω, the population likelihood L can be calculated by the individual likelihood L by Equation 1e. i is associated with.

number

number

number

[0031] Common probabilistic programming languages ​​such as Stan or Turing can be used to determine these posterior distributions using Markov Chain Monte Carlo (MCMC) techniques such as Hamiltonian Monte Carlo. Estimated distribution plots of NLME parameters calculated using Pumas and other software via HMC on the data show how discrete MCMC chains can be reinterpreted as probability distributions via kernel density estimation (KDE).

[0032] 1.2 Method Once the model parameters and their uncertainties are determined, the model is used to determine an individualized dose schedule. The dose schedule specifies how much, where, and when to administer a drug to a subject. As used herein, "dose" refers to the amount of drug administered to a specific location within a subject at a specific time, not the resulting concentration level in the subject's target tissue. Due to a random component that accounts for inter-individual variability and event-to-event variability of results within a single individual, a particular dose will yield a range of possible outcomes in terms of the drug concentration actually observed. Understanding the range of possible outcomes should inform the decision regarding the original dose. In many previous approaches, the expected outcome is expressed in terms of the probability of an actual concentration occurring within the therapeutic range. If that probability is high enough for a particular dose, then that dose is considered acceptable.

[0033] It should be noted that a subset of all acceptable doses may be preferred because the probability of harmful deviation from the therapeutic range is smaller than other administered doses within the acceptable range. Thus, in various embodiments, a cost function is defined that includes a measure of weighted deviation from the therapeutic range, and the dose determined is based on reducing or minimizing the weighted deviation. This formulation reduces or minimizes the cost, which is based on the probability of high deviation of the subject's drug response from a known safety window. Although this cost function may yield high values, it is a better indicator of the actual cost to the subject and is therefore more useful for optimizing the dose (compared to minimizing the distance from a particular optimal response curve, which does not capture the fact that a solution within the safety window is "good"). Thus, the quantity of interest in some embodiments is the probability of rare but large deviations.

[0034] In some embodiments, this cost is evaluated based on determining the distribution of initial conditions, e.g., doses, using the Koopman adjoint approach, which is several times more efficient than a naive Monte Carlo approach for determining that distribution. In some embodiments, the increased efficiency is achieved by employing, in addition to the Koopman adjoint approach, data-parallel numerical solutions to ordinary differential equations (ODEs) that are accelerated on a graphics processing unit (GPU).

[0035] Bayesian estimation or other techniques such as prior knowledge or maximum likelihood can be used to estimate θ and η i It is assumed that a sufficient estimate of the posterior distribution of x is obtained to answer the following question: What is the optimal or advantageous dosing regimen Di for patient i? The remainder of this section focuses on this single patient i, and therefore the subscript i is dropped in this section.

[0036] To illustrate a typical approach, we introduce the Frobenius-Perron (FP) operator. Define S(x)=u(t), where u(0)=x. In other words, S(x) is the solution of a dynamical system with initial conditions x and t is time. The FP operator is used to compute a probabilistic statement of the patient outcome with respect to the uncertainty in the initial conditions. First, note that the uncertainty in the parameters (θ and η) is incorporated into this calculation. For ordinary differential equations, the expanded system of ODEs is expressed in Equation 1b (repeated below as Equation 2), Equation 3, and Equation 4.

number

number

number

number

number

number

[0037] Given a dynamical system

number

number

[0038] For dose, the population uncertainty is fixed and the individual random effects, i.e., θ and η with probability distribution f(x), induce a probability distribution in the solution of the dynamics u(t), and pushing the probability distribution through this dynamical system is expressed as the FP operator, i.e.,

number

[0039] Regarding the issue of dosage, see P SIt may be useful to know the inverse action of u = u ≠ u ≠ y ...

[0040] Koopman operator U S is defined by Equation 6.

number

number

number

[0041] 2B is a set of plots illustrating the interaction of uncertainty in initial conditions, represented by pdfs, with a cost function using a Koopman adjoint process, according to one embodiment. In plots 230 and 260, the horizontal axis shows the initial conditions x, and the vertical axis is dimensionless, ranging from 0 to 1. The assumed known pdf of the initial conditions f(x) is illustrated by dashed trace 236. In plots 240 and 250, the horizontal axis shows the resulting measurement variable y, and the vertical axis is dimensionless, also ranging from 0 to 1, representing the probability density of y. In this space of plots 240 and 250, the cost function g(y) is illustrated by solid trace 247. The FP operator transforms the pdf of the initial conditions trace 236 into a result pdf, illustrated by dashed trace 246 in plot 240. The dot product of these two, representing the left hand side of Equation 8, is given by the shaded region enclosed by trace 248. This is the result pdf weighted by the cost function. Expected value E[g(X)|X~P S f] is represented by the area of ​​the shaded region. S is a function of the cost function g, illustrated by trace 247 of plot 250, over the initial condition space U, illustrated by trace 267 of plot 260. S g form. U, as illustrated by trace 236, which represents the right-hand side of Eq. S The dot product of g and the initial condition pdff is given by the shaded region enclosed by trace 268. This is the pdf of the initial conditions weighted by the cost function. The expectation E[U S g(X)|X~f] is represented by the area of ​​the shaded region. The dot products represented by the areas of the filled regions under traces 248 and 269 are equal, respectively.

[0042] Given the uncertainty f(x) in the propagation of initial conditions x through the dynamical system S, the optimal dose is the dose that has the highest expectation of a good patient outcome. In other words, we wish to find a dose schedule D that optimizes the expectation with respect to a cost function g at a solution u (or measurement y) of the dynamical system S. Let g(u(t)) be the cost associated with a given patient outcome, then this is the preferred dose plan D that satisfies Equation 9. * corresponds to finding

number

number

number

[0043] The advantage of using Equation 11 can be understood as follows: A common way to perform dose optimization with respect to the uncertainty of the parameters represented in Equation 9 is to use P S The solution to this problem is to use Monte Carlo estimation of f, which involves 1) finding sample parameters θ and ηi from the uncertainty distribution f, 2) solving the dynamical system to compute S(x) at each set of parameters, and 3) computing g(S(x)) at each set of parameters and taking the average individually. This procedure is computationally expensive since it requires the solution of many differential equations. Using Equation 11, the desired expectation can be calculated much more efficiently, for example using multidimensional quadrature, where U Sg(x) is the solution of the dynamical system with parameters determined by a quadrature procedure, and f(x) is the evaluation of the probability distribution at the quadrature points. The quadrature procedure is explained in the document Cubature_(Multi-dimensional integration), in subfolder index.php of folder wiki in subdomain initio of domain mit of super-domain edu.

[0044] The cost function g is advantageously based on clinical trials used to establish safety guidelines known as the therapeutic window. For example, for a given drug, the area under the curve (total concentration or AUC) of the drug concentration in the target tissue or other dynamic variable u (or a measure y that depends on u) is known to be safe over 24 hours when 200-400 milligrams is the typical amount at which the therapeutic window is written. Given how every individual metabolizes at a given rate, the goal is to optimize the probability that the dose is within the therapeutic window with respect to the uncertainty of the patient-specific effect. Any cost function g(u) can be used. In some embodiments, the cost function not only accounts for "goodness" in giving preferential weight (lower cost) to results within the therapeutic window, but also allows for the cost of "badness," such as avoiding certain ranges by assigning higher costs to certain ranges or increasing costs for distance from the therapeutic window. The use of distance outside the therapeutic window encompasses the idea that being close to the therapeutic window satisfies the condition, while being farther away may pose greater risk to the patient.

[0045] To illustrate this method, one embodiment using a simple cost function is described herein: g(y)=χw(y), where χw is an indicator function of the therapeutic window W, e.g., g is a function defined as 1 if a measurement y of a solution u to a differential equation with parameters (θ, η, D, Z) falls within the therapeutic window W, and 0 otherwise. In this definition of the cost function, the expected value of the characteristic function is the probability that an event occurs, and thus this definition of g leads to Equation 12.

number

[0046] Numerical methods for computing integral expressions such as Equation 11 are called quadrature methods. For the computation of the quadrature, the method can use several efficient computational techniques. In embodiments where the dimension of the integral (the number of elements in the initial vector x) is small enough, techniques such as Cubature can be used effectively. However, in other embodiments where the dimensionality of the parameter space is higher, dimension-independent quadrature methods such as Monte Carlo importance sampling methods like VEGAS can be used to increase the computational efficiency. The naive Monte Carlo solution is given by the central limit theorem as

number

[0047] This reformulation can significantly reduce the computational cost of computing the cost functions, since the computation of each integrand requires only one solution of the dynamical system, which is the most expensive part of the computation (usually the numerical solution to the ODE). Furthermore, each of these quadrature methods allows for a high degree of parallelism. Implementations such as the Julia programs Cuba.jl, available from the folder cuba at domain feynarts of super-domain de, and Cubature.jl, available from the folder ￣stevenj at subdomain math of domain mit of super-domain edu, allow for simultaneous integrand computations of a user-selected batch size M.

[0048] In one embodiment, a Julia language tool such as DiffEqGPU.jl can be called to automatically generate a data-parallel GPU-parallelized ODE solver from any existing pure Julia ODE solver code with generic AbstractArray processing using a non-standard interpretation. The ODE integrator generated by this meta-program simultaneously computes the solution of ODEs of various parameters, splitting the computation among many cores of the GPU, and the solution is then used to compute the loss at all M parameters. As this is true for any existing pure Julia ODE solver, the method generates GPU-based integrators for both non-stiff and stiff ODEs, including implicit-explicit (IMEX) methods and adaptive methods for stochastic differential equations (SDEs), delay differential equations (DDEs), differential algebraic equations (DAEs), hybrid differential equations such as jump-diffusion or Markov regime switching models, totaling over 200 methods. The loss for each of the M trajectories is then weighted by the probability of the parameters to complete the computation of the integrand.

[0049] 3 is a flow diagram illustrating an exemplary method for determining a dose that will minimize cost, according to one embodiment. Although the steps are shown in FIG. 3 as integral steps in a particular order for illustrative purposes, in other embodiments, one or more steps or portions thereof may be performed in a different order or overlapping in time, in series or in parallel, or omitted, or one or more additional steps may be added, or the method may be modified in any number of combinations.

[0050] In step 301, a NLME model is developed using appropriate random and non-random parameters, and model parameters are learned based on, for example, historical data from an EHR. Training continues until a certain stopping condition is reached, such as a number of iterations, or the maximum or percentile difference between the model output and the measurements is less than a certain tolerance, or the difference in results between successive iterations is less than a certain target tolerance, or the results diverge, or some other known measure, or some combination. As a result of training, the model is a continuous multivariate model that does not need to perform classification or binning of subjects as required in the application of prior art nomograms. The model includes at least one parameter based on a population mean, at least one parameter (e.g., Ω) that represents random deviations between individuals, and at least one part that represents dynamics (e.g., using an ordinary differential equation that represents changes over time). As a result of step 301, a continuous multivariate model is received at the processor as first data indicating at least one distribution parameter for dose-response to a pharmaceutical drug.

[0051] In step 311, a cost function is determined based on a safety range, such as a therapeutic range or a therapeutic index or some combination. In some embodiments, the cost function includes a weighting function of distance from the therapeutic range. For example, results within the therapeutic range have no cost, and increasing distance beyond the therapeutic range increases the cost. This can be linear or have any shape function depending on the distance. In some embodiments, particularly toxic ranges, for example as determined based on the therapeutic index, are weighted higher, not necessarily proportional to the distance from the therapeutic range. In an embodiment implemented automatically by the processor, step 311 includes receiving second data indicative of a therapeutic range for the dose response, and receiving third data indicative of a cost function based on the distance of the dose response from the therapeutic range.

[0052] In step 313, the Koopman transform of the cost function is determined or programmed on the processor, for example by defining a numerical procedure to calculate g(S(x)), which defines a function that receives the initial condition x0 and the parameters of the population including the distribution parameters, perform the dynamics, and evaluate the cost.

[0053] In step 315, for the candidate doses, an expected cost of the candidate dose is determined based on at least a parameter of random deviation (e.g., Ω for each of one or more model parameters) and a Koopman transformation of the cost function using multi-dimensional quadrature. In some of these embodiments, the evaluation is performed based on, for example, known graphics processing unit acceleration techniques.

[0054] In step 321, it is determined whether the cost is below a particular target threshold cost, such as a cost found to be associated with a good outcome. If not, control proceeds to step 323 to adjust the candidate dose, for example, by decreasing the candidate dose by an incremental amount or percentage. Control then returns to step 315 to recalculate the cost. If the assessed cost is below the threshold, control proceeds to step 331.

[0055] At step 331, the subject is administered the drug for each candidate dose. Thus, when the expected cost is below the threshold cost, the candidate dose of the drug will be administered to the subject.

[0056] In step 333, samples are taken from the subject to determine the therapeutic effect (e.g., concentration of drug in the sampled tissue) at measurement times equal to or greater than zero determined during step 325. This step is performed in a timely manner to correct shortcomings in the administration of the drug before adverse effects, such as ineffective or toxic outcomes, occur in the subject. This step contrasts with previous methods in which a subsequent dose of drug is administered regardless of whether the first dose was ineffective or toxic to the individual subject. Thus, step 327 involves sampling the subject in one or more sets to obtain corresponding values ​​of corresponding measures of therapeutic effect.

[0057] In step 335, the model determined in step 305 is updated based on the measurements made in step 333, if any. For example, if the drug concentration value is less than or greater than the therapeutic range, the individual is determined to have a response that deviates substantially from the population mean. Any method can be used to update the model parameters to fit the observed response of the individual, such as new training of a neural network. In some embodiments, the population model is updated based on the observation of one of the individuals being treated in proportion to the number of different subjects that previously contributed to the model. In some embodiments, in addition to or instead of updating the population model parameter values, model parameter values ​​specific to a particular subject are generated.

[0058] Thus, method 300 determines an individualized protocol for administering medications that are likely to remain within therapeutic range in a subject during treatment of the subject. Additional formats are useful because they provide significant performance and feature advantages.

[0059] As personalized precision dosing becomes more commonplace, efficient and robust methods for optimizing doses with respect to probabilistic outcomes will become more central to the practice of pharmacometrics. As personalized computing becomes more prevalent, more accurate models will be needed, and computing will need to be transferred to patients living with embedded or mobile devices that have little computing power.

[0060] 2. Exemplary Embodiments In various exemplary embodiments, the dose is determined using some or all steps of method 300 for the pharmaceuticals listed in the table below. For each pharmaceutical listed, the concentration window W or therapeutic index is listed in the third column. The therapeutic index (TI, also called therapeutic ratio) is a quantitative measure of the relative safety of a drug. It is a comparison of the amount of therapeutic drug that produces a therapeutic effect to the amount that produces toxicity. The cost function g(y) proposed herein is based, at least in part, on the values ​​in this third column. Koopman expectation is used in Equation 11 or Equation 12 with a quadrature procedure to estimate the optimal range or subset of the dose ranges listed in the second column. [Table 1] [Table 2] [Table 3] [Table 4] [Table 5] [Table 6] [Table 7]

[0061] According to an exemplary embodiment, at least a portion of the method 300 is performed to determine and administer a dosing regimen to deliver and maintain blood levels of theophylline within the therapeutic range. This embodiment demonstrates how the use of Koopman Expectation can accelerate the calculations 100x over Monte Carlo while providing numerical error bounds to improve clarity of the results. The methodology is implemented using the high performance differential equation solver DifferentialEquations.jl in the Pumas Pharmacometrics Suite, so that existing models can be automatically adapted with the accelerated calculations. Being a pure Julia software stack, the methodology can be compiled into mobile devices such as ARM for direct deployment to patients.

[0062] This exemplary embodiment is a high-performance implementation of Koopman expectation in Pumas. Pumas is a high-performance pharmaceutical modeling and simulation engine written in the Julia language. At its core, Pumas uses the DifferentialEquations.jl ecosystem to drive event-based models that can be used to characterize the time course of drugs or other endogenous substances in the human body. Such characterization in Pumas facilitates accurate dose prediction and optimization of treatment trajectories. The program was developed in Julia programming (in the domain julialang super-domain org). Equation 11 is calculated using DifferentialEquations.jl [3] as U given a Pumas specification of a dynamical system S. Sg was implemented by computing the . In other embodiments, alternative ODE solvers for dynamic equations such as Sundials (available at subdomain computing of domain llnl of super-domain gov in folder projects subfolder sundails.gov subfolder projects document sundials) or LSODA (available at subdomain people of subdomain sc of subdomain fsu of super-domain edu in folder ￣jburkardt subfolder f77_src subfolder odepack document odepack.html) can also be used to evaluate the dynamics. Multidimensional integrals were calculated using Quadrature.jl, a wrapper library for popular quadrature methods such as Cuba and Cubature mentioned above. This integral implementation allows for a batch solver that parallelizes the calculations at the quadrature points, allowing for quadrature to be multithreaded, distributed, and GPU accelerated. Table 8 provides code that demonstrates the use of Koopman expectation to calculate the probability of an AUC less than 300 for the theophylline model. The definition of the nonlinear mixed-effects model is defined using the standard @model macro format in Pumas. In this example, there are no individual random effects η, and the @param block defines the fixed effect θ as theta, indicating that it is a vector of four variables. The second blog, @covariates, defines the covariates as Z i =[sex i ,wt i ,etn i ], where the specific value for subject i is obtained from the patient data. The pre block calculates the dynamic system [Ka i ,CL i ,V i In this example, we determine the value of the random effect η i Note that the is omitted. @dynmics is the dynamic equation

number

[0063] Next, we load a data set of values ​​for the covariate Z to define the subjects and dosing schedule. Then, we calculate θ and η iWe input the distribution of the uncertainty of the parameters of . This would normally arise from a Bayesian posterior probability distribution as described above, but for illustrative purposes a prescribed distribution is specified. For demonstration purposes we use θ~[Uniform(1.0,5.0),Uniform(20.0,100.0),Uniform(200.0,400.0),Uniform(0.1,2.0)], where θ denotes a four-dimensional random variable whose first value is defined by a uniform random variable of Uniform(1.0,5.0), and similarly for the other values. The probability distribution of Uniform(a,b) is f(x)=1 / (ba) for x∈[a,b] and zero otherwise. We then define the therapeutic region in terms of the loss function g of the observables. For this implementation, we use Pumas NLME integration and NCA to calculate the AUC as part of the variables derived from the model, and therefore the AUC exists as one of the observed outputs of the system. We can therefore define the therapeutic region as below a threshold. Here we call the expectation routine and tell it to use the Koopman Expectation method. We specify that we use the HCubeJL quadrature method, which tends to be efficient for low-dimensional integrand functions (<8 uncertainty variables). This calculation then calculates a given expectation value up to a desired tolerance. Given its differentiability, this function can be used in an optimization loop with a defined dose plan, thus performing dose optimization under uncertainty.

[0064] To demonstrate the significant performance gains and feature benefits, the KoopmanExpectation option was compared to the MonteCarloExpectation option, and the convergence speed of the two methods was determined for calculating the probability of a patient outcome occurring within the therapeutic window, with various imaxiser choices (number of ODE solver allowances) increased. FIG. 4 is a graph showing an example convergence of a Koopman Expectation calculation compared to a previous method using Monte Carlo Expectation, according to an example embodiment. The horizontal axis shows the number of ODE solver calls, which is a measure of computational cost, and the vertical axis shows the probability estimate. Here, Koopman Expectation was evaluated using routine Cubature with the integration shown as trace 406. The convergence of Monte Carlo is shown as trace 408. Both are given in terms of the number of ODE solver calls required. Note that the cubature integration method trace 408 is accompanied by a free error estimate, shown by bounds 407.

[0065] Figure 4 demonstrates two results. One result is that the Koopman method converges to provide stable probability estimates in the last two orders of magnitude with nearly 1 in 100 ODE solutions, and the computational cost is roughly two orders of magnitude less. The second result is that by using the H-Cubature method, the Koopman calculation not only determines the probability, but also gives numerical error bounds on the probability estimates, allowing the user to effectively know the uncertainty introduced by numerical errors. Such bounds are very difficult to compute with Monte Carlo estimation, given the slow convergence rate of the variance. Taken together, this demonstrates the Koopman expectation as both an efficient and robust method.

[0066] 3. Overview of Computational Hardware FIG. 5 is a block diagram illustrating a computer system 500 on which an embodiment of the present invention can be implemented. The computer system 500 includes a communication mechanism, such as a bus 510, for passing information between other internal and external components of the computer system 500. The information is represented as a measurable phenomenon, usually a physical signal of voltage, but in other embodiments includes phenomena such as magnetic, electromagnetic, pressure, chemical, molecular atomic, and quantum interactions. For example, a north and south magnetic field, or zero and non-zero voltages represent the two states (0, 1) of a binary digit (bit). Other phenomena can represent numbers of a higher base. A superposition of multiple simultaneous quantum states before measurement represents a quantum bit (qubit). A series of one or more digits constitutes digital data, which is used to represent numbers or codes of letters. In some embodiments, information, referred to as analog data, is represented by an approximate continuum of measurable values ​​within a certain range. The computer system 500, or portions thereof, constitute a means for performing one or more steps of one or more methods described herein.

[0067] The series of binary digits constitutes digital data used to represent numbers or codes for characters. Bus 510 contains many parallel conductors of information so that information may be transferred quickly between devices coupled to bus 510. One or more processors 502 are coupled to bus 510 for processing information. Processor 502 performs a set of operations on the information. The set of operations includes taking information from bus 510 and putting information onto bus 510. The set of operations also typically includes comparing two or more units of information, shifting the position of the information units, and combining two or more units of information, such as by addition or multiplication. The series of operations performed by processor 502 constitutes computer instructions.

[0068] Computer system 500 also includes a memory 504 coupled to bus 510. Memory 504, such as random access memory (RAM) or other dynamic storage device, stores information including computer instructions. Dynamic memory allows the information stored therein to be changed by computer system 500. RAM allows units of information stored at locations called memory addresses to be stored and retrieved independently of information at adjacent addresses. Memory 504 is also used by processor 502 to store temporary values ​​during execution of computer instructions. Computer system 500 also includes read only memory (ROM) 506 or other static storage device coupled to bus 510 for storing static information including instructions that are not changed by computer system 500. Also coupled to bus 510 is a non-volatile (permanent) storage device 508, such as a magnetic or optical disk, for storing information including instructions that persists even when computer system 500 is turned off or loses power.

[0069] Information, including instructions, is provided to the bus 510 for use by the processor from external input devices 512, such as a keyboard with alphanumeric keys operated by a human user or sensors. The sensors detect conditions in their vicinity and convert those detections into signals compatible with those used to represent information within the computer system 500. Other external devices coupled to the bus 510, used primarily for human interaction, include a display device 514, such as a cathode ray tube (CRT) or liquid crystal display (LCD), for presenting images, and a pointing device 516, such as a mouse or trackball or cursor direction keys, for controlling the position of a small cursor image presented on the display 514 and for issuing commands related to graphical elements presented on the display 514.

[0070] In the illustrated embodiment, special purpose hardware such as application specific integrated circuits (ICs) 520 are coupled to bus 510. The special purpose hardware is configured to perform operations not performed by processor 502 fast enough for a special purpose. Examples of application specific ICs include graphics accelerator cards for generating images for display 514, cryptography boards for encrypting and decrypting messages sent over a network, voice recognition, and interfaces to specialized external devices such as robotic arms and medical scanning equipment that repeatedly perform certain complex sequences of actions that are more efficiently implemented in hardware.

[0071] The computer system 500 also includes one or more instances of a communication interface 570 coupled to the bus 510. The communication interface 570 provides a two-way communication coupling to various external devices, such as printers, scanners, and external disks, that operate with their own processors. Typically, the coupling is to a network link 578 that connects to a local network 580 to which various external devices with their own processors are connected. For example, the communication interface 570 may be a parallel port or a serial port or a universal serial bus (USB) port on a personal computer. In some embodiments, the communication interface 570 is an integrated services digital network (ISDN) card or a digital subscriber line (DSL) card or a telephone modem that provides a telecommunications connection to a corresponding type of telephone line. In some embodiments, the communication interface 570 is a cable modem that converts signals on the bus 510 to signals for a telecommunications connection via a coaxial cable or optical signals for a telecommunications connection via a fiber optic cable. As another example, the communication interface 570 may be a local area network (LAN) card for providing a data communication connection to a compatible LAN, such as an Ethernet. Wireless links may also be implemented. Carrier waves, such as sound and electromagnetic waves, including radio, optical, and infrared, travel through space without wires or cables. Signals include man-made variations in the amplitude, frequency, phase, polarization, or other physical characteristics of the carrier wave. In the case of a wireless link, communications interface 570 transmits and receives electrical, acoustic, or electromagnetic signals, including infrared and optical signals, that carry information streams, such as digital data.

[0072] The term computer-readable medium is used herein to refer to any medium that participates in providing information to the processor 502, including instructions for execution. Such media may take many forms, including but not limited to non-volatile media, volatile media, and transmission media. Non-volatile media include, for example, optical or magnetic disks, such as the storage device 508. Volatile media include, for example, dynamic memory 504. Transmission media include, for example, coaxial cables, copper wire, fiber optic cables, and waves that travel through space without wires or cables, such as acoustic and electromagnetic waves, including radio, optical, and infrared. The term computer-readable storage medium is used herein to refer to any medium that participates in providing information to the processor 502, except for transmission media.

[0073] Common forms of computer readable media include, for example, a floppy disk, a flexible disk, a hard disk, a magnetic tape or any other magnetic medium, a compact disk ROM (CD-ROM), a digital video disk (DVD) or any other optical medium, a punch card, a paper tape or any other physical medium having a pattern of holes, a RAM, a programmable ROM (PROM), an erasable PROM (EPROM), a FLASH-EPROM, or any other memory chip or cartridge, a carrier wave, or any other medium from which a computer can read. The term non-transitory computer readable storage medium is used herein to refer to any medium that participates in providing information to the processor 502, excluding carrier waves and other signals.

[0074] Logic encoded in one or more tangible media includes one or both of computer-readable storage media and processor instructions on dedicated hardware such as ASIC 520.

[0075] Network link 578 typically provides information communication through one or more networks to other devices that use or process the information. For example, network link 578 may provide a connection through a local network 580 to a host computer 582 or to equipment 584 operated by an Internet Service Provider (ISP). The ISP equipment 584 in turn provides data communication services through the public, world-wide packet-switched communication network of networks now commonly referred to as the Internet 590. A computer connected to the Internet, referred to as a server 592, provides services in response to information received over the Internet. For example, server 592 provides information representing video data for presentation on display 514.

[0076] The present invention relates to the use of computer system 500 to implement the techniques described herein. According to one embodiment of the present invention, these techniques are performed by computer system 500 in response to processor 502 executing one or more sequences of one or more instructions stored in memory 504. Such instructions, also referred to as software and program code, may be loaded into memory 504 from another computer-readable medium, such as storage device 508. Execution of the sequences of instructions stored in memory 504 causes processor 502 to perform the method steps described herein. In alternative embodiments, hardware, such as application specific integrated circuit 520, may be used in place of or in combination with software to implement the present invention. Thus, embodiments of the present invention are not limited to any specific combination of hardware and software.

[0077] Signals transmitted through communication interface 570 over network link 578 and other networks carry information to and from computer system 500. Computer system 500 can transmit and receive information, including program code, through networks 580, 590, among others, through network link 578 and communication interface 570. In an example using Internet 590, server 592 transmits program code for a particular application requested by a message transmitted from computer 500 through Internet 590, ISP equipment 584, local network 580, and communication interface 570. The received code may be executed by processor 502 as received, or may be stored in memory device 508 or other non-volatile storage for later execution, or both. In this manner, computer system 500 can obtain application program code in the form of signals on a carrier wave.

[0078] Various forms of computer readable media may be involved in carrying one or more sequences of instructions and / or data to the processor 502 for execution. For example, the instructions and data may initially be carried on a magnetic disk of a remote computer, such as the host 582. The remote computer loads the instructions and data into its dynamic memory and sends the instructions and data over a telephone line using a modem. A modem local to the computer system 500 receives the instructions and data on the telephone line and uses an infrared transmitter to convert the instructions and data into a carrier wave infrared signal that acts as the network link 578. An infrared detector, acting as the communications interface 570, receives the instructions and data carried in the infrared signal and places information representing the instructions and data on the bus 510. The bus 510 carries the information to the memory 504, from which the processor 502 retrieves and executes the instructions, using some of the data transmitted with the instructions. The instructions and data received by the memory 504 may optionally be stored on storage device 508 either before or after execution by the processor 502.

[0079] FIG. 6 illustrates a chipset 600 in which an embodiment of the present invention can be implemented. The chipset 600 includes a processor and memory components, as described in connection with FIG. 5, programmed to perform one or more steps of the methods described herein, e.g., integrated into one or more physical packages (e.g., chips). By way of example, the physical package includes an arrangement of one or more materials, components, and / or wires on a structural assembly (e.g., baseboard) to provide one or more characteristics, such as physical strength, maintaining size, and / or limiting electrical interactions. In certain embodiments, it is contemplated that the chipset can be implemented on a single chip. The chipset 600, or portions thereof, constitute a means for performing one or more steps of the methods described herein.

[0080] In one embodiment, chipset 600 includes a communication mechanism, such as a bus 601, for passing information between components of chipset 600. Processor 603 has connectivity to bus 601 for executing instructions and processing information stored, for example, in memory 605. Processor 603 can include one or more processing cores, each configured to execute independently. A multi-core processor allows for multi-processing within a single physical package. Examples of multi-core processors include two, four, eight, or more processing cores. Alternatively or additionally, processor 603 can include one or more microprocessors configured in tandem via bus 601 to allow for independent execution of instructions, pipelining, and multi-threading. Processor 603 can also include one or more specialized components for performing specific processing functions and tasks, such as one or more digital signal processors (DSPs) 607 or one or more application specific integrated circuits (ASICs) 609. DSPs 607 are typically configured to process real-world signals (e.g., sound) in real time, independent of processor 603. Similarly, ASIC 609 may be configured to perform specialized functions not easily performed by a general purpose processor. Other specialized components that assist in carrying out the functions of the invention described herein may include one or more field programmable gate arrays (FPGAs) (not shown), one or more controllers (not shown), or one or more other special purpose computer chips.

[0081] The processor 603 and associated components are connected to memory 605 via bus 601. Memory 605 includes both dynamic memory (e.g., RAM, magnetic disks, writable optical disks, etc.) and static memory (e.g., ROM, CD-ROM, etc.) for storing executable instructions that, when executed, perform one or more steps of the methods described herein. Memory 605 also stores data associated with or generated by the execution of one or more steps of the methods described herein.

[0082] 4. Substitutions, Deviations and Modifications In the above specification, the present invention has been described with reference to specific embodiments thereof. However, it will be apparent that various modifications and changes can be made thereto without departing from the broader spirit and scope of the present invention. The specification and drawings are therefore to be regarded in an illustrative and not a restrictive sense. Throughout this specification and claims, unless the context requires otherwise, the word "comprise" and its variations such as "comprises" and "comprising" are understood to imply the inclusion of the stated item, element, or step or group of items, elements, or steps, but not the exclusion of any other item, element, or step or group of items, elements, or steps. Furthermore, the indefinite article "a" or "an" is intended to indicate one or more of the items, elements, or steps modified by the article.

[0083] 5.References The following citations are referenced above: [Table 9]

Claims

1. 1. A method automatically executed by a processor for generating an individualized dosage protocol, comprising: receiving first data indicative of a continuous multivariate model of a population with respect to dose response to a pharmaceutical agent, the model having at least one distribution parameter characterizing the variability of the population; receiving second data indicative of a therapeutic range for said dose response; receiving third data indicative of a cost function based on a distance of the dose-response from the therapeutic value range; assessing expected costs for candidate doses based at least in part on the distribution parameters and a Koopman transformation of the cost function; administering the candidate dose of pharmaceutical agent to a subject when the expected cost is below a threshold cost; The method includes:

2. The method of claim 1 , wherein evaluating the expected costs further comprises performing a multi-dimensional quadrature.

3. The method of claim 1 , wherein the at least one distribution parameter is based on a distribution of values ​​in an electronic health record (EHR).

4. 10. The method of claim 1, wherein the pharmaceutical agent is of a type selected from the group including antiepileptic drugs, anticoagulants, antiarrhythmic drugs, and immunosuppressants.

5. The method of claim 1 , wherein the distance of a dose-response from the therapeutic value range of the cost function is indicated by a therapeutic index.

6. 1. A non-transitory computer readable medium carrying one or more instruction sequences for generating an individualized dosage protocol, the one or more instruction sequences upon execution by one or more processors causing the one or more processors to: receiving first data indicative of a continuous multivariate model of a population with respect to dose response to a pharmaceutical agent, the model having at least one distribution parameter characterizing the variability of the population; receiving second data indicative of a therapeutic range for said dose response; receiving third data indicative of a cost function based on the distance of the dose-response from the therapeutic value range; assessing expected costs for candidate doses based at least in part on the distribution parameters and a Koopman transformation of the cost function; administering the candidate dose of the pharmaceutical agent to a subject; A non-transitory computer-readable medium having a program for executing the program.

7. The computer-readable medium of claim 6 , wherein evaluating the expected costs further comprises performing a multi-dimensional quadrature.

8. The computer-readable medium of claim 6 , wherein the at least one distribution parameter is based on a distribution of values ​​in an electronic health record (EHR).

9. 7. The computer readable medium of claim 6, wherein the pharmaceutical agent is of a type selected from the group including antiepileptic drugs, anticoagulants, antiarrhythmic drugs, and immunosuppressants.

10. 7. The computer-readable medium of claim 6, wherein the distance of the dose-response from the therapeutic value range of the cost function is indicated by a therapeutic index.

11. 1. An apparatus for generating an individualized dosage protocol, comprising: At least one processor; at least one memory containing one or more sequences of instructions; wherein the at least one memory and the one or more instruction sequences are provided to the at least one processor to provide the device with at least receiving first data indicative of a continuous multivariate model of a population with respect to dose response to a pharmaceutical agent, the model having at least one distribution parameter characterizing the variability of the population; receiving second data indicative of a therapeutic range for said dose response; receiving third data indicative of a cost function based on a distance of the dose-response from the therapeutic value range; assessing expected costs for candidate doses based at least in part on the distribution parameters and a Koopman transformation of the cost function; administering the candidate dose of pharmaceutical agent to a subject when the expected cost is below a threshold cost; An apparatus configured to cause the

12. The apparatus of claim 11 , wherein evaluating the expected costs further comprises performing multi-dimensional quadrature.

13. The apparatus of claim 11 , wherein the at least one distribution parameter is based on a distribution of values ​​in an electronic health record (EHR).

14. 12. The apparatus of claim 11, wherein the medicine is of a type selected from the group including antiepileptic drugs, anticoagulants, antiarrhythmic drugs, and immunosuppressants.

15. The apparatus of claim 11 , wherein the distance of a dose-response from the therapeutic value range of the cost function is indicated by a therapeutic index.

Images