From: https://towardsdatascience.com/predicting-pharmacokinetics-with-deterministic-models-and-bayesian-statistics-e3822a28977b

Introduction

Pharmacokinetics

Pharmacokinetics (PK) deals with the distribution and metabolism of xenobioticsin organisms. In drug development and clinical research, one encounters pharmacologists whose job it is to advise physicians the dose of drug(s) to provide to their patients. For this purpose, they often run simulations of so-called pharmacokinetic models which are deterministic, i.e. they consist of ordinary differential equationswith respect to a time variable and that describe the following processes (known as LADME):

• liberation: the process of drug release from the pharmaceutical formulation
• absorption: the process of a substance entering the blood circulation
• distribution: the dispersion or dissemination of substances throughout the fluids and tissues of the body
• metabolism: the recognition by the organism that a foreign substance is present and the irreversible transformation of parent compounds into daughter metabolites
• elimination: the removal of the substances from the body

Such models help determining the optimal dose of drugs as a function of patient attributes like weight, height, age, or sex. In addition, the concentration profile of a drug in the blood plasma over time can be predicted. Furthermore, drug interactions can be monitored. This way, the dosing interval can be defined for an individual patient — paving the way for personalized treatments(visit openPK for more information).

Deterministic Modeling

Normally, these sets of equations that are part of deterministic models have parameters which are unknown and have to be estimated from experimental data. The body itself is described as a multi-compartment system (organs are compartments) in which subsets of the mentioned processes take place, and drug/metabolite molecules flow from one compartment (organ) to another. It is usually assumed that said compartments are ideally mixed, i.e. incoming molecules are immediately mixed. Note that time-dependent experimental data is usually not available for all compartments, especially not in humans.

Model Development

Model Formulation

The way to develop such a model is to formulate a set of differential equations (one example dealing with ethanol metabolism programmed in MATLAB can be found on my GitHub). The equation below is a general formulation of a mass balance with a flow rate Q and a volume V for a compound i in a compartment k. P is a partition coefficient which states how much a compound is retained in a certain tissue. Q and V depend on patient weight, height, age, and sex (modeled with allometric scaling).

dc_dt_i_k = (Q_k/V_k)*(c_in_i_(k-1)-c_i_k/P_i_k)

Values for Q, V,and P for each compartment can be found in the literature. In addition, drugs are metabolized in the liver compartment and excreted in the kidney compartment. In this case, the equations for the liver and kidneys are modified.

dc_dt_i_liver = (1/V_liver)*(Q_liver*(c_in_i-c_i_liver/P_i_liver)-v_max_i*c_i_liver/(K_m_i+c_i_liver))dc_dt_i_kidney = (1/V_kidney)*(Q_kidney*(c_in_i-c_i_kidney/P_i_kidney)-k_r_i*c_i_kidney

The term in bold in the liver equation represents the so-called Michaelis-Menten expression for a reaction catalyzed by an enzyme. It limits a reaction to a maximal reaction speed vmax since the amount of active sites for reaction is limited. However, vmax and Km are not readily available in the literature and have to be estimated from experimental data (the same is true for the clearance term kr). Furthermore, competition for active sites with other molecules might occur and the reaction term might need a modification by including an inhibition expression. It is a priori not clear which processes to include in which compartments. For example, a drug is absorbed into the blood within the gastrointestinal tract, but it could also decompose there. The way to go is to formulate an initial model, to estimate its parameters, and to evalute its performance on a test set. If the model performance is insufficient, the model is refined by adding or removing some terms, and the procedure is repeated until convergence.

Parameter Estimation

There are different objective functions to be used for the parameter estimation. First, the regular least-squares problemcan be minimized (the concentrations are a function of parameters p and so is the objective function J). The concentrations of all compounds i in all compartments k are stacked in one column for this purpose (using a reshape function to avoid a for-loop). The variable c_m refers to the measured concentration, whereas c_p is the predicted concentration.

J = (1/N)*(c_m-c_p)'*(c_m-c_p)

This function has the disadvantage that it will fit molecules in high concentration better than such of low concentrations as their contribution to the objective is much smaller. Another possibility is to scale the dot product with the inverse of a covariance matrix and reformulate it as a maximum likelihood problem. Of course, the variance matrix E has to be known.

J = (1/N)*(c_m-c_p)'*inv(E)*(c_m-c_p)

E can, for example, be estimated from the coefficient of variation cv for each compound in each compartment.

E = cv'*c_m

Finally, by defining a prior distribution over the parameters, the objective function can be expanded and the problem solved in a maximum a posteriori manner, which is a method from Bayesian statistics. In this case, the objective function is again slightly modified assuming normal distribution with zero mean and qI (q times identity matrix variance and neglecting constants).

J = (1/N)*(c_m-c_p)'*inv(E)*(c_m-c_p)+q*p'p

Lastly, the model has to be wrapped inside the objective function. With the parameters, the model is evaluated and the predicted concentrations are calculated using a solver for differential equations (e.g. ode15s to deal with stiffness), i.e. by a Runge-Kutta method. The objective function is minimized using a solver like fminunc which is based on a Quasi-Newton methodimplementation (similar to gradient descent, but the gradient is approximated around a current parameter point and scaled proportionally to the magnitude of the second derivate of the objective function, which enhances the performance by allowing bigger steps in flat regions). It’s important to perform a multi-start (different starting points for parameters p) in order to find a global optimum and not to end up with a local one. (Note that a least-squares solution can serve as a starting point for a more sophisticated approach.)

Model Validation

The model is best validated using a test set, i.e. data with which the model has not been trained. Mean relative/absolute (squared) error can be used as metrics. In addition, it is important to perform a sensitivity analysis. For this purpose, the partial derivatives of all concentrations (every compound in every compartment) with respect to all parameters are calculated. Note that they are time-dependent as the concentrations are time-dependent. Hence, the influence of one parameter on an output variable can change in time. Very sensitive parameters and the processes they describe have to be further studied.

The rank of the sensitivity matrix reveals if all parameters are independent of each other (structural identifiability); if the rank is lower than the number of parameters, the model has to be redefined by removing or readjusting dependent parameters.

R = rank(S)

Additionally, the confidence intervals of the parameters should be calculated to examine practical identifiability. The variance of the parameters can be calculated using the Fisher informationtogether with the Cramér-Rao bound.

E_p = inv(S'*inv(E)*S)

The confidence intervals are calculated by a univariate approach using only the diagonal elements of the parameter covariance matrix and the quantiles of Student’s t-distribution.

sigma = sqrt(diag(E_p))CI_lb = p_fit-sigma*tinv(0.95, number_points-number_parameters)CI_ub = p_fit+sigma*tinv(0.95, number_points-number_parameters)

If a confidence interval includes zero, then a parameter could also be zero, i.e. not taking part in the equations, hence the model should be refined. Note that all parameters should be positive.

After several iterations of fitting and refinement, the construction process is certainly finished, a reasonable test set performance is obtained, and the model can be used for prediction.