Input Estimation for Extended-Release Formulations Exemplified with Exenatide

Estimating the in vivo absorption profile of a drug is essential when developing extended-release medications. Such estimates can be obtained by measuring plasma concentrations over time and inferring the absorption from a model of the drug’s pharmacokinetics. Of particular interest is to predict the bioavailability—the fraction of the drug that is absorbed and enters the systemic circulation. This paper presents a framework for addressing this class of estimation problems and gives advice on the choice of method. In parametric methods, a model is constructed for the absorption process, which can be difficult when the absorption has a complicated profile. Here, we place emphasis on non-parametric methods that avoid making strong assumptions about the absorption. A modern estimation method that can address very general input-estimation problems has previously been presented. In this method, the absorption profile is modeled as a stochastic process, which is estimated using Markov chain Monte Carlo techniques. The applicability of this method for extended-release formulation development is evaluated by analyzing a dataset of Bydureon, an injectable extended-release suspension formulation of exenatide, a GLP-1 receptor agonist for treating diabetes. This drug is known to have non-linear pharmacokinetics. Its plasma concentration profile exhibits multiple peaks, something that can make parametric modeling challenging, but poses no major difficulties for non-parametric methods. The method is also validated on synthetic data, exploring the effects of sampling and noise on the accuracy of the estimates.


STRUCTURAL IDENTIFIABILITY ANALYSIS
The structural identifiability of two PK models has been considered, Gao and Jusko (2012) and Li et al. (2015).

Structural identifiability analysis approaches
The structural identifiability analysis in this paper has been done using two approaches: The Exact Arithmetic Rank (EAR) and the Taylor series expansion approach.
The EAR approach was published in Karlsson et al. (2012) and for details of the method the reader is referred to that publication. The method is implemented as an automated Mathematica package and considers structural local identifiability. The user defines the model states, model parameters, initial conditions and input and output functions. These are used as arguments to the structural identifiability function defined within the Mathematica package which returns whether the model is locally identifiable or not. The method itself is based on local algebraic observability (Sedoglavic, 2002). Higher order derivatives of the output function(s) are generated resulting in a system of algebraic equations. The inverse function theorem is applied by checking whether the Jacobian matrix of the algebraic system is rank deficient or otherwise. If the Jacobian matrix of the algebraic system has full rank the model is structurally identifiable. All computations are carried out numerically, meaning that the model parameters and the initial conditions are numerically instantiated with randomly chosen integers. The input function is approximated with a truncated power series expansion with integer coefficients.
The Taylor series approach was first published in Pohjanpalo (1978) and for details of the method the reader is referred to that publication. This approach makes use of the fact that the coefficients in a power series expansion of a function uniquely determines the trajectory of that function. By computing the power series expansion of the output function and determine whether the model parameters can be uniquely determined from the coefficients in the power series expansion the structural identifiability of the model can be studied.

Gao and Jusko model
The structural identifiability analysis was done only of the PK model in Gao and Jusko (2012). In equation (1) in the paper, the variable R is used. This was interpreted as the total amount of free receptors at time t, i.e., receptors available for drug binding at time t. Therefore, in the structural identifiability analysis the variable R was replaced with (R tot − RC(t)).
The model structure of the PK model is given by where the unknown parameter vector is {V c , k el , k pt , k tp , k on , k of f , k int , R tot }, C(t) is the drug concentration in the central compartment, A T (t) is the drug amount in a peripheral compartment, RC(t) is the concentration of the drug-receptor complex, and u(t) is the unknown input.

Structural identifiability results
These were the results following from the structural identifiability analysis using the EAR approach.
• The PK model with bolus IV dose is structurally identifiable. • The PK model with IV infusion is structurally identifiable.
• The PK model with SC administration is structurally unidentifiable with one degree of freedom. The unidentifiable parameters are F and V c , the rest of the model parameters are identifiable.
Since all of the model parameters are shared for the IV case to the SC case, measures that can be taken to make the SC model identifiable include: i) fixing the parameter V c to the estimated value from the IV case.
ii) setting up a joint inference problem with both models.
In the paper the bioavailability parameter F was fixed to a value which solves the structural identifiability issue. However, this was not done as a direct result from a structural identifiability analysis.

Identifiability of the input function
In the analysis of whether the input function is identifiable all model parameters were assumed to be known.
To analyze whether the input function in the PK model in Gao and Jusko (2012) is identifiable or not the Taylor series expansion approach was used (Pohjanpalo, 1978). A Taylor series expansion of a function y(t) around t = 0 is given by where k = 1, 2, . . . All coefficients in the Taylor series determines uniquely the function y(t). By computing the Taylor series expansion around t = 0 of the output function we obtain the following coefficients: The unknowns in the equations above are {u(0),u(0),ü(0), u (3) (0), . . . } since the model parameters are assumed to be known. Because of the structure of this particular model the u (n−1) (0) term will always enter linearly in the expression (see last term of each equation) for y (n) (0). As can be seen above, this results in a triangular structure of the higher order derivatives and we can therefore conclude that {u(0),u(0),ü(0), u (3) (0), . . . } can be determined from the output function.
The Taylor series expansion of the input signal u(t) around t = 0 is where k = 1, 2, . . . Since we have shown that all of the coefficients in the Taylor series expansion of u(t) can be determined we can conclude that the input function u(t) is identifiable.

SYNTHETIC INPUT FUNCTION
The chosen synthetic input function u(t) is a sum of three functions: where u 1 (t) and u 2 (t) are given by an Erlang distribution: and u 3 (t) is given by a biexponential function: The functions u 1 (t) and u 2 (t) can each be interpreted as the output of a chain of n i compartments, connected with a rate constant of k tr i , where the initial value of the first compartment is a i and all other compartments are initialized to 0. This rewriting is known as the "Linear Chain Trick" (Jacquez, 1996;Smith, 2010). Similarly, the initial biexponential u 3 (t) can be interpreted as the output from an absorption compartment combined with a peripheral compartment. The complete input function can be generated by the model in Fig. S1. The parameters k a , k 12 and k 21 in the figure are related to r 1 , r 2 , k 1 and k 2 via: where a 0 is the initial amount in compartment 1, and q is an intermediate variable to simplify the expressions. This result can be obtained by analytically solving the system of differential equations corresponding to compartments 1 and 2, and matching the coefficients. It might be more biologically plausible to connect u 1 (t) and u 2 (t) to the absorption compartment rather than directly to the central compartment. However, this would make little difference, since the parameter k a is an order of magnitude larger than k tr 1 and k tr 2 , while making the algebraic expressions more complicated. Figure S1. Compartmental model interpretation of the functions generating the test data that were used for assessing the performance of the estimation method. Compartments 1-2 generate the fast biexponential function accounting for the initial peak. The two chains 3-12 and 13-16 each produce an Erlang distributionshaped function accounting for the peaks seen at longer timescales. Double circles indicate compartments with non-zero initial conditions. C is the central compartment. Note that compartments 6-10 are not shown.

ABSORPTION RATE PLOTS
The paper shows the estimation results for the real exenatide data expressed as the absorbed fraction of the drug as a function of time, as this is how such results are typically presented in the context of ER-formulation development. This appendix shows the same estimation results expressed as the absorption rate as a function of time, which can be useful for identifying the phases of drug absorption. 10.0 mg exenatide ER Plasma concentration Figure S2. Estimation results including 95% credible intervals for the initial 48 hours of the real exenatide data, shown as absorption rate. It can be seen that the variations in the data for the 7.0 mg dose causes a larger uncertainty region. 10.0 mg exenatide ER Plasma concentration Figure S3. Estimation results including 95% credible intervals for the whole 12-week period of the real exenatide data, shown as absorption rate.