
==== Front
J Pharmacokinet Pharmacodyn
J Pharmacokinet Pharmacodyn
Journal of Pharmacokinetics and Pharmacodynamics
1567-567X
1573-8744
Springer US New York

37864654
9887
10.1007/s10928-023-09887-3
Original Paper
Learning pharmacometric covariate model structures with symbolic regression networks
Wahlquist Ylva ylva.wahlquist@control.lth.se

Sundell Jesper
Soltesz Kristian
https://ror.org/012a77v79 grid.4514.4 0000 0001 0930 2361 Department of Automatic Control, Lund University, P.O. Box 118, 221 00 Lund, Sweden
21 10 2023
21 10 2023
2024
51 2 155167
14 6 2023
18 9 2023
© The Author(s) 2023
2023
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Efficiently finding covariate model structures that minimize the need for random effects to describe pharmacological data is challenging. The standard approach focuses on identification of relevant covariates, and present methodology lacks tools for automatic identification of covariate model structures. Although neural networks could potentially be used to approximate covariate-parameter relationships, such approximations are not human-readable and come at the risk of poor generalizability due to high model complexity.In the present study, a novel methodology for the simultaneous selection of covariate model structure and optimization of its parameters is proposed. It is based on symbolic regression, posed as an optimization problem with a smooth loss function. This enables training of the model through back-propagation using efficient gradient computations.Feasibility and effectiveness are demonstrated by application to a clinical pharmacokinetic data set for propofol, containing infusion and blood sample time series from 1031 individuals. The resulting model is compared to a published state-of-the-art model for the same data set. Our methodology finds a covariate model structure and corresponding parameter values with a slightly better fit, while relying on notably fewer covariates than the state-of-the-art model. Unlike contemporary practice, finding the covariate model structure is achieved without an iterative procedure involving manual interactions.

Keywords

Pharmacometrics
Covariate modeling
Pharmacokinetics
Symbolic regression
Neural networks
Knut and Alice Wallenberg FoundationELLIIT Strategic Research AreaLund UniversityOpen access funding provided by Lund University.

issue-copyright-statement© Springer Science+Business Media, LLC, part of Springer Nature 2024
==== Body
pmcIntroduction

Pharmacokinetics (PK) are the dynamics governing drug uptake, distribution and elimination from the body. For many drugs, it is common to assert a low-order linear and time-invariant (LTI) compartment model for PK modeling. Such low-order models typically capture the uptake, distribution, and elimination dynamics adequately. Increasing model complexity through, for example, additional compartments results in models where the parameters are not practically identifiable from data collected during clinical trials or clinical practice.

While suitable PK model structures (e.g. number of compartments and topology) for common drugs can be established from data, a notable challenge exists in that the parameter values that explain the PK for one individual, are often not suitable for another. Such inter-individual variability is partly explainable by individual-specific features that are referred to as covariates, and partly attributed to random effects.

The purpose of pharmacometric covariate modeling is to identify covariates (i.e., fixed effects) responsible for inter-individual variability, thus minimizing random effects. The gold standard is to approach this problem in a Bayesian setting using mixed-effect modeling, for which there exists both mature [1] and novel [2] tools. However, to apply mixed-effect modeling, one needs to decide which of possibly many covariates (e.g. age, body mass, gender, or genetic factors) to include. One also needs to decide on a parametric function that maps included covariates to parameters of the PK model. There is a tradition of using functions that are of sufficiently low complexity to be human-readable, and some function classes are more popular than others [3, 4]. However, for data sets including numerous covariates, selection of covariates and functions to consider is a combinatorial problem, where the search space may become limiting. Furthermore, due to step-wise approaches typically used for covariate identification, functions including multiple covariates are typically only considered if supported by prior knowledge [5].

Recently, the interest in combining machine learning (ML) and pharmacometrics [6] has increased. For example, ML has shown promising results in concentration predictions [7], identification of influential covariates [8, 9], as well as parameter regression and model selection with the use of genetic algorithms and neural networks (NNs) [10]. The ML methods have been able to match, and even beat, classic NLME modeling at a much higher computational speed. However, using ML to select influential covariates and identify the structural model simultaneously has, to our knowledge, not yet been studied.

In the present paper, we provide an automatic method for simultaneous covariate selection and identification of covariate functions using an adaptation of ML. Similarly to the current standard approach, we maintain simple human-readable expressions. The method is based on symbolic regression [11], that approximates the underlying combinatorial problem with a smooth continuous one, thus enabling the use of efficient gradient-based optimization methods.

To illustrate utility, we apply our methodology to a large data set for the drug propofol and compare the resulting covariate model with the current state-of-the-art model [4]. Our method produces a model with increased predictive performance, using fewer covariates. Furthermore, the covariate model is optimized autonomously in contrast to the prevalent iterative and manual procedure.

Methods

The method described in this paper aims to automatically learn closed-form pharmacometric parameter–covariate relationships from data. We will consider learning expressions for PK model parameters from drug administration and resulting blood plasma concentration data, both of which are time series, but not necessarily synchronous or equidistant samples. The overarching setup for this is shown in Fig. 1 and we will focus on a concrete example based on a large multi-study data set published in [4] for the anesthetic drug propofol, commonly modeled with a mammillary three-compartment model [12].Fig. 1 Mammillary three-compartment model example illustrating our novel method. The objective is to automatically learn the covariate model f that maps a known covariate vector φ (comprising e.g., age, gender, or genetic factors) to the parameter vector θ (e.g. rate constant k·· and volumes V·) of a fixed-structure pharmacometric (PK) model. The method is data-driven in that it uses drug administration profiles (time series data) u, and model-based in that it assumes a PK model of known structure. In this example, f is learned to minimize some error measure between observed (i.e. measured from samples) blood plasma concentrations and corresponding predictions Cpred by the model. Dots in the graphs show instances of dose changes and blood samples, respectively

Data set

We demonstrate our method on a data set composed of propofol plasma concentration observations from 1031 individuals1 from 30 clinical studies aggregated by Eleveld et al. [4], from here on referred to as the data set. Ethical approvals of the underlying studies are declared in the original publications, referenced by [4]. The data set contains 15,433 observations, of which 11,530 were arterial and 3903 venous. Of the 1031 individuals, there were 670 males and 361 females with ages ranging from 27 weeks to 88 years, and weights ranging from 0.68 to 160 kg, see Table 1.Table 1 Covariate candidates considered in the PK model of [4] as well as in our modeling. Individuals of the data set fall within the reported ranges

Covariate	Interpretation	Unit	Range	
φ1	Age	Years	0–88	
φ2	Weight	kg	0.68–160	
φ3	BMI	kg m0-2	6.2–52.8	
φ4	Gender	Male/female	670 M/ 361 F	
φ5	Blood sampling site	Arterial/venous	727 A/ 306 V	

The main reason for choosing this data set as a demonstrator is that propofol is a drug with well-studied pharmacokinetics. Evidence of this, which also constitutes a benchmark for our demonstrator, is the model presented in [4]. Furthermore, the data set is that it has been openly disclosed by Eleveld et al., enabling transparent third-party analysis of our work.

In our model development, we consider age, weight, BMI, gender, and blood sampling site (arterial or venous) as potential covariates of our PK model. These are the same candidates as considered in [4], where individual demographics were disclosed as part of the data set.

The data was pre-processed the same way as in [4]: data points corresponding to subsequent infusion changes spaced closer than 1 s apart in the time dimension, or 0.5 μgs-1 in the dose dimension were merged.

Pharmacokinetic model

We consider a three-compartment mammillary model to describe the pharmacokinetics of propofol. The drug concentration xi [μgL-1] in compartment i∈{1,2,3} is 1a x˙1=-(k10+k12+k13)x1+k21x2+k31x3+1V1u,

1b x˙2=k12x1-k21x2,

1c x˙3=k13x1-k31x3,

where kij describes the drug transfer rate [1/s] from compartment i to j. The drug is administered at rate u [μgs-1] to the central compartment (i=1), which is also where the propofol plasma concentration is measured. The volume of the central compartment is V1 [L].

In the literature, the equivalent parameterization of volumes (V1,V2,V3) and clearances (CL,Q2,Q3) constitutes a common alternative to Eq. 1. The conversions between these parameterizations are 2a CL=k10V1[Ls-1],

2b Q2=k12V1[Ls-1],

2c Q3=k13V1[Ls-1],

2d V2=k12k21V1[L],

2e V3=k13k31V1[L].

We have chosen to implement our method using the parameterization in Eq. 1 due to numeric benefits. These are further explained in [13], where we developed a fast and natively differentiable simulator for the three-order mammillary model.

Covariate model

The covariate model, shown in Fig. 1, is expressed as a function f that maps the covariate vector φ=[φ1,…,φnφ]⊤ to the vector θ=[θ1,…,θnθ]⊤ of PK model parameters. Thus f has components f1,⋯,fnθ, each mapping the covariate vector φ to one of the nθ PK model parameters.

In our example with the three-compartment model for propofol, there are nφ=5 covariates and nθ=6 PK model parameters according to Table 1. To not favor covariates based on their scale, all input covariates were normalized during model development. Continuous inputs (age, weight and BMI) were scaled from 0 to 1, and categorical inputs (gender and blood sampling site) were scaled to ±0.5.

Predictive performance

To assess the quality of a particular covariate model candidate f, we need a performance measure that captures how well f reflects the training data. Most optimization methods, including the ones used in this paper, relies on this measure being scalar.

For comparability, we employ the same scalar performance measures as those used in [4]: We train our model to minimize an ensemble average of absolute logarithmic error (ALE). Subsequently, we evaluate predictive performance in terms of ALE, and three additional error measures: logarithmic error (LE), prediction error (PE), and absolute prediction error (APE).

For each individual in the data set, there is a vector of observation–prediction errors, where each entry corresponds to the difference between a blood sample observation and the value predicted by the model. Observation-prediction errors could thus be computed for each sample, over a time series for one individual, or the entire data set. This prompts a consistent notation and we index by ij the error over a sample j for individual i, whereas a total error over a time series for an individual is indexed by i.

The per-sample (absolute) logarithmic error (A)LE [14] is thus 3a LEij=ln(Cobsij/Cpredij),

3b ALEij=|LEij|,

where Cobsij are observed (measured) plasma concentrations and Cpredij are corresponding predictions produced by the model.

Similarly, the per-sample (absolute) prediction error (A)PE [15] is 4a PEij=Cobsij-CpredijCpredij·100%

4b APEij=|PEij|.

To avoid divisions by zero if Cpredij=0, a special casing is needed where the error is set to zero for such samples.

Using Eqs. 3 and 4, corresponding per-individual errors were computed by taking the median, resulting in the median absolute logarithmic error (MdALE):5 MdALEi=medianALEij,j=1,⋯,ni,

where ni is the number of entries of the time series from individual i.

To train the model of [4] and the ones considered here, a single scalar loss function representing model fit is required. This is obtained by averaging across the individuals. Training the covariate model to minimize ALE thus translates into minimizing6 JALE=1n∑i=1nMdALEi,

where n is the number of individuals in the data set.

Similarly to the analysis performed in [4], ALE and APE were used as indicators of model accuracy and LE and PE were used as indicators of bias. Values closer to zero for ALE and APE reflect better accuracy and values closer to zero for LE and PE indicate less bias.

Clinically acceptable ranges for MdPE0i and MdAPE0i are 10-20% and 20-40%, respectively [4, 15, 16]. This translates to acceptable clinical ranges for bias measure MdLE0i<0.18 and accuracy measure MdALE0i<0.34.

Of note, the choice of loss function will affect how outliers are penalized, where using an average loss across individuals would for example penalize outliers more than a median loss over individuals.

Symbolic regression networks

At the core of our methodology lies a small artificial neural network (ANN) with a specific structure, that of a symbolic regression network [17], representing a simple closed-form expression. In our case, the purpose is to learn human-readable closed-form expressions describing the covariate model f. A schematic illustration of such network is shown in Fig. 2.Fig. 2 Symbolic regression network with three layers, each marked by a gray box. The output of node zli at layer l is the ith component of zl=Wlxl+bl, where xl, Wl and bl are the input vector, weight matrix, and bias vector of that layer. The base expressions gli acting on zli take on the role of activation functions used in ordinary ANNs. For example, the output of the first layer (and therefore input to the second layer) is x2=g1(W1x1+b1). Input and output of the network is the covariate vector φ=x1 and the PK parameter θk=x4, respectively

In general, ANNs can be viewed as flexible function approximators that can be trained to represent input–output mappings that fit available data. An ANN consists of nl layers, where the output vector xl+1 of layer l is obtained by applying a vector of nonlinear activation functions gl to an affine transformation zl to the layer input xl: 7a zl=Wlxl+bl,

7b xl+1=gl(zl).

The linear weight matrices Wl and bias vectors bl constitute the free parameters used to train the network. In the conventional case, the components of the activation functions gl are monotonously increasing functions, such as sigmoids [18]. In contrast, a symbolic regression network can be understood as an ANN with the additional constraint that the resulting approximator f should be a human-readable closed-form expression of a mathematical function [17, 19–21]. Training of symbolic regression networks is therefore often referred to as equation learning [17]. The methodology has gained broad attention for its ability to produce impressive results in discovering (known) physical laws from data [22].

The activation functions g of a symbolic regression network represent mathematical functions that we refer to as base expressions. In standard symbolic regression, a sequence of topologies with different base expressions and scalar parameters–of which the weight matrices and bias vectors in Eq. 7 would constitute a special case–are evaluated in search of one that fits the available data well [23]. Conventional equation learning is thus a combinatorial problem, with poor (exponential) time complexity, and therefore often approached using genetic algorithms [23].

Instead of relying on random perturbations as in genetic algorithms, we use gradient-based optimization methods common to conventional ANN training. To ensure human-readability, we enforce sparsity of the ANN by alternating training epochs with pruning epochs, in which the least important parameters are removed from the network. We next rely on a concrete example based on the Eleveld data set, to describe and illustrate this approach.

We set out with a nominal (unpruned) network of Fig. 2. It constitutes our nominal representation of fk, where the inputs are the covariates according to Table 1. Thus, the output of the network represents one of the PK parameters θk, such as k10 or V1 of Eq. 1. For a model with nθ PK parameters, we use a parallel interconnection of nθ symbolic regression networks, each modeling one component fk of f, mapping the covariate vector φ to each PK parameter θk,k=1,⋯,nθ.

In our example, the base expressions of each layer l∈{1,2,3}, with input vector zl=zl1,zl2,⋯⊤, were chosen to cover previously published PK models for propofol, such as [4] and [12]:g1(z1)=z11z12·z13|z14|z15,g2(z2)=z21z22·z23z24z25+1,g3(z3)=|z3|,

where the g3 assures positive output of the final layer. The division in g2 has the term one in the denominator to assure that the output does not blow if z25≥0 approaches zero.

Training

Minimizing the loss function Eq. 6 across trainable parameters of the covariate model f is referred to as training. In our case, these parameters are stacked into a vector γ, made up by the elements of all weight matrices and bias vectors. Since the data set is static, we can view training as minimization of the scalar-valued function JALE(γ). In order to do so, we need to evaluate JALE(γ). Doing so requires simulation of one PK model for each data set individual, to obtain the predicted plasma concentration at each observation time instance. This can be efficiently done using the method introduced in [13].

In addition to the supporting fast simulation, the method of [13] enables exact and efficient evaluation of derivatives of individual prediction errors with respect to the trainable model parameters. This allows us to train the covariate model f using conventional ANN back-propagation, using the stochastic gradient-based optimization algorithm ADAM [24]. As with artificial neural networks in general, establishing formal convergence guarantees is challenging. However, in practice both standard deep learning networks and our symbolic regression do converge to models that fit data adequately well. Similarly to the deep learning case, convergence rate will also vary with data, choice of activation functions, learning rate, and initialization of the training. Our implementation, using the neural network package Flux [25], relies on the differential programming capabilities of the Julia language [26]. A full disclosure of our implementation is found in the GitHub repository [27].

Pruning

To obtain simple human-readable expressions from an initially dense symbolic regression network, such as the one in Fig. 2, we alternate between parameter training of the fixed network structure and pruning of the network.

The goal of this pruning is to obtain a sparse network structure, which translates into a readable expression of the corresponding covariate model. In the pruning process, we remove covariates and network parameters that have relatively little influence on the network output. This process is visually exemplified in in Fig. 3.

It is desirable to only include as many covariates as needed to explain the data [5]. Therefore, we start by identifying and removing the least important covariates from the symbolic regression network to obtain expressions with fewer covariates. Next, we prune network parameters γ (linear weights and biases) to obtain simple covariate expressions. Pruning a network parameter is achieved by fixing its value to zero and removing it from the vector γ of trainable parameters.

Before each pruning iteration, we train the network until convergence. This translates into finding a local minimum of the loss function in Eq. 6. At such minima, the partial derivatives of the loss function with respect to the trainable parameters are zero. Therefore, the second derivatives define the first non-zero terms of the Taylor series expansion that describe local parameter sensitivities, as further explained in Appendix A. The second-order derivatives make up the elements of the Hessian matrix. Specifically its diagonal elements represent sensitivities in the corresponding individual parameters.

In the Taylor series expansion, the second-order terms take on the form8 S(γk)=γk2Hk,

where Hk is the kth Hessian diagonal element of the loss function with respect to the network parameter γk. S(γk) denotes salience of a parameter γk [28].

For pruning of the network parameters γ, Eq. 8 may be used directly. However, computing the salience for a covariate is not as straightforward since the Hessian elements would differ between individuals. A solution to this is to sum the salience contribution of each individual,9 S(φk)=∑i=1nφik2(Hi)k,

where (Hi)k is the kth diagonal element of the Hessian, Hi, determining the sensitivity of the covariate φik for individual i, and n is the number of individuals.

We use the Zygote package [29] in Julia to compute the Hessian diagonal elements. In Appendix A, we give a more detailed description of the Hessian-based pruning method, providing mathematical insight into this methodology.Fig. 3 Pruning sequence of a symbolic regression network with output pharmacokinetic parameter θ2=k12. Input covariates are age, weight, gender, and arterial or venous sampling (AV). The nominal network has three dense layers, each followed by base expressions such as 1 (feedforward), multiplication, power function, division and absolute value. A black line represents a connection between two nodes, and a gray line represents a pruned (removed) connection. The final network represents the covariate expression of Eq. 10b

Recipe: Symbolic regression

A compact summary of our training and pruning scheme is provided below. Initially, we start with a nominal symbolic regression network, like the one shown in Fig. 2. It is sequentially trained and pruned until we obtain a final expression of sufficient complexity and fit. Our training and pruning sequence of a symbolic regression network is as follows: Choose a nominal symbolic regression network architecture, and corresponding base expressions.

Train the network until convergence.

Compute the salience of each (remaining) covariate, S(φk) of Eq. 9.

Sort the covariates by saliency and remove the covariate with the smallest salience.

If the desired final number of covariates is reached, continue to step 6, otherwise return to step 2.

Train the reduced network until convergence.

Compute the salience of each (remaining) trainable network parameter, S(γk) of Eq. 8.

Sort the network parameters by saliency and remove N parameters with the smallest salience.

If the desired final number of network parameters is reached, continue to step 10, otherwise return to step 6.

Train the reduced network until convergence.

Convert the resulting network to a readable functional expression.

The number of covariates and network parameters to keep in the final symbolic regression network is a trade-off between fit to data and complexity of the final expression. In our example, illustrated in Fig. 3, we remove one covariate at each pruning iteration, until only two covariates are left. Next, we remove the N=10 least sensitive parameters in the first parameter pruning iteration, and then one parameter per subsequent pruning iteration until only twelve network parameters are left. More details of the pruning and training can be found in [27]. In [30], we demonstrate that our methodology can identify functions of known shapes correctly.

Initialization of the parameter values associated with step 1 of the recipe affects the fit of the resulting final model. To mitigate the risk of poor model fit due to an unfortunate initialization, the recipe could be run several times, where the best fitting model is kept. In our example, we have executed the recipe eight times.

Limits of performance

Structural mismatch between the asserted PK model structure Eq. 1 and the data, in combination with measurement errors, induce upper and lower limits on prediction errors of the trained covariate model. Here we explain how these limits can be characterized by training two additional models.

Even with a very complex covariate model, one cannot expect perfect fit to data (zero loss). This is because the fixed (three-compartment) PK model structure is only an approximation of the actual pharmacokinetics, combined with (blood sample) measurement errors. Part of the loss remaining after applying our covariate modeling scheme can thus be attributed to this mismatch. To indicate how much, we optimize one set of (three-compartment) PK parameters for each individual in the data set. This results in a completely covariate-free model. While it will fit data better than any covariate model, it does not generalize. This makes it practically useless for purposes other than providing an upper limit for performance.

Another natural question to ask is how much we gain (in terms of loss) by considering covariate dependencies. To do this, we optimize a constant, i.e. covariate-free, model—one where all individuals share the same (three-compartment) PK parameter values. This model does constitute a lower bound of the achievable predictive performance that any covariate-based model should beat.

For these two additional models, PK parameter optimization was done with the optimization package Optim.jl [31], to minimize the same loss as for our covariate model based on symbolic regression. Implementation details can be found in [27].

Results

Applying the proposed methodology to the Eleveld data set [4], resulted in the following covariate model, mapping covariates to rate constants and central compartment volume of the PK model in Eq. 1: 10a k10=0.00441WGTWGTmax+0.00342[s-1]

10b k12=0.158AGEAGEmax-0.00431BMIBMImax-0.1880.64AGEAGEmax-0.0174BMIBMImax-0.743[s-1]

10c k13,male=0.0058AGEAGEmax2+0.00208AGEAGEmax+0.00262.75AGEAGEmax2+0.985AGEAGEmax+0.601[s-1]

10d k13,female=0.0058AGEAGEmax2-0.00208AGEAGEmax+0.00262.75AGEAGEmax2-0.985AGEAGEmax+0.601[s-1]

10e k21=0.00408BMIBMImax2-8.16·10-4BMIBMImax-0.0057BMIBMImaxWGTWGTmax+0.00218[s-1]

10f k31,male=4.52·10-5+1.92·10-5AGEAGEmax[s-1]

10g k31,female=4.52·10-5-1.92·10-5AGEAGEmax[s-1]

10h V1=0.0596AGEAGEmax+18.7WGTWGTmax-13.7WGTWGTmax2-3.5AGEAGEmaxWGTWGTmax-0.0557[L].

AGE, WGT, BMI represent age [years], weight [kg] and body mass index [kg m-1]. The subscript max represent the input normalization where AGEmax=88 years, WGTmax=160 kg and BMImax=52.8 kg m-1.

The subscript male or female indicates different PK parameter expressions depending on gender. The blood sampling site (arterial or venous) was available as a modeling covariate, but was automatically pruned by the symbolic regression algorithm.

The obtained model of Eq. 10 is less complex than the Eleveld model in [4], provided in Appendix B for reference, and comparable to simpler covariate models for propofol, such as [3, 16].

The predicted concentrations of the final covariate model on the Eleveld data set are shown in Fig. 4 together with the corresponding Eleveld model predictions. The distribution of errors between predicted and observed predictions is shown in the boxplot in Fig. 5, indicating comparable predictive capability of the models, despite our model being less complex, and involving fewer of the available covariates. For ease of comparison, we present the corresponding average prediction errors in Table 2. We compare our model to the Eleveld model, to a covariate-free PK model where all individuals share the same parameters values, and to covariate-free individual models with unique parameter sets. As expected, adding covariates explains some, but not all, variability between patients. This can be seen by comparing our model to the constant PK model and the individual models. As seen in Table 2, all of the prediction errors: MdLE, MdALE, MdPE, and MdAPE fall within clinically acceptable ranges, as presented further above.Fig. 4 Predicted versus observed propofol concentrations of our covariate model (Symreg, red) compared to the Eleveld covariate model in [4] (Eleveld, blue) in logarithmic scale. The identity function, representing a perfect model fit, is shown in black

Fig. 5 Comparison of prediction error MdALE Eq. 5 between predicted and observed propofol concentrations for pharmacokinetic models. Our covariate model is denoted Symreg and the Eleveld covariate model is described in [4]. The constant model represents one parameter set over the population and individual represents individual set of model parameters. The lower whisker for the individual models goes to zero

Table 2 Comparison of prediction errors for the propofol data set in [4] (1,031 individuals) for several pharmacokinetic models, all trained with MdALE Eq. 5 as loss. Individual and constant PK model(s) are shown for comparison, representing best and worst case limits, respectively

Method	Mean MdALE	Mean MdLE	Mean MdAPE	Mean MdPE	
Constant model	0.501	0.140	211	181	
Eleveld model	0.325	0.0791	34.8	14.4	
Symbolic regression	0.279	-0.0489	27.4	-0.266	
Individual models	0.0623	1.65 ·10-4	6.21	0.0178	

Discussion

We have introduced a novel symbolic regression methodology for simultaneously automating the search for a suitable covariate model structure and optimization of its parameters. Similar to contemporary methodologies, it relies on a user-specified set of base expressions from which the covariant model can be composed. However, and in important contrast to contemporary methods, the need for a combinatorial search across combinations of these expressions is voided.

Throughout the paper, a propofol PK modeling example has been used as a demonstrator. Within this example, the proposed methodology manages to accomplish slightly better data fit than the result of state-of-the-art modeling [4] (see Figs. 4 and 5, Table 2), while relying on fewer covariates, see Eq. 10. The obtained values of individual volumes and clearances were comparable to those in [4].

The introduced methodology is broadly applicable to PK modeling from time series data, and likewise to pharmacodynamic (PD) modeling, and combined pharmacokinetic and pharmacodynamic (PK/PD) modeling. It’s main benefit lies in that it poses the search for a suitable PK model as symbolic regression with a smooth loss function. This, in combination with efficient methods for simulation and gradient computations [13] enables efficient model learning using back-propagation.

Another advantage is that the method can find covariate functions that are both simple and explain available data, while having a structure that would generally not be considered in a manual model structure search, unless explicitly supported by prior knowledge. However, the obtained covariate model is deterministic in the sense that we do not obtain a distribution over individual PK parameter values. It would be possible to integrate this methodology into the traditional mixed-effect modeling framework. Yet, this would be computationally much more demanding, which is why we propose first applying symbolic regression in a deterministic setting to arrive at a covariate model structure, and then (if desired) apply mixed-effect modeling to maximize parameter likelihood (with respect to some parameter priors and subject to the considered data) within the found structure.

In this paper, we have focused on models that are of sufficiently low complexity to be human-readable, which is achieved by enforcing sparsity of the neural network that constitutes the expression tree of the covariate model. If human readability is not necessary, an ordinary deep ANN could be employed instead. However, there is also a trade-off between fit to training data (expressiveness) and generalization to yet unseen data due to possible over-fitting. Enforcing human-readability naturally limits flexibility of the model, thus decreasing this risk of over-fitting. The ability to manually specify base expression also provides a means to integrate expert knowledge into the model. For example, a suspicion that a compartment volume should correlate to the square of patient age would motivate multiplication as a base expression. This enables incorporating expert knowledge into the model, for example a clearance to the power of 0.75 as in [32], or compartmental allometry as in [33].

We assessed the model’s out-of-sample prediction performance using a five-fold cross-validation. The data set was divided into five equal parts, and a covariate model was trained on four of the partitions, excluding one each time. The excluded partition, referred to as the validation set, was used to evaluate the model. This process was repeated for all five partitions to obtain an average predictive performance. The resulting mean (range) MdALE from cross-validation on the test sets was 0.303 (0.152--0.423), similar to the mean MdALE of 0.279 in table 2. There was thus only a modest 10 % difference in mean error between the training and validation sets, which suggests no over-fitting issue. If such a problem had occurred, reducing the number of covariates and parameters (harder pruning) could have balanced the prediction errors on both sets.

The way we have employed the methodology here differs from how covariate models are usually trained using mixed-effect modeling. Rather than asserting parameter priors and selecting the most likely parameters from the posterior distribution that the data infers, we have chosen to train a scalar loss function, resulting in a deterministic covariate model. It would in theory be possible to embed our methodology within an inference engine, but at the cost of high computational cost. If a Bayesian interpretation is desired, a likely better alternative is to first run our methodology to arrive at a covariate model structure–as we have done in our example–and then assert parameter priors to the parameters of the resulting model, to finally apply mixed-effect modeling to compute the corresponding posteriors.

Conclusion

We have presented a novel methodology for automatic and simultaneous covariate model structure discovery and parameter optimization. This model was demonstrated using an example on which it outperforms state-of-the-art modeling, in that it finds expressions that match data slightly better, while relying on notably fewer covariates. We conclude that the potential of automated model structure discovery is substantial; it could greatly optimize the process of pharmacometric covariate modeling. Additionally, it’s likely to provide an improved balance between model complexity and data fit. This improved balance is something that could be challenging to achieve when simply assessing a series of pre-set model structure candidates in sequence.

Appendix A. Hessian-based pruning

In this appendix we provide a mathematical interpretation of the roles played by the Hessian and parameter saliences in our pruning approach.

The effect of perturbing the parameter vector can be analyzed by approximating the loss function J(γ) by a Taylor series. A perturbation ∂γ of the parameter vector γ will change the loss function by∂J(γ)=∇J(γ)∂γ+12∑iHii∂γi2+12∑i≠jHij∂γi∂γj+O(||∂γ||3),

where ∂γi are the components of ∂γ, ∇J(γ) is the gradient of the loss function and Hij are the elements of the Hessian matrix H of J with respect to γ so that∇J(γ)=∂J(γ)∂γHij=∂2J(γ)∂γi∂γj.

The goal is to prune the parameters that are least sensitive, i.e. those that affect the loss J the least when perturbed. Even for our size of network, repeatedly computing the full Hessian H would notably slow down training of the symbolic regression model. Instead, we make a diagonal approximation of the Hessian, neglecting cross terms Hij with i≠j.

Allowing training to converge before each pruning iteration, ensures that ∇J(γ)=0 (or in practice negligible), thus enabling approximation of ∂J(γ) by∂J(γ)≈12∑iHii∂γi2.

This approximation is then used to form an importance measure (salience) of our network parameters, according to [28], where the salience of parameter γi becomesS(γi)=Hiiγi2.

Appendix B. The Eleveld model

The Eleveld PK covariate model in [4] is given here,fageing(x)=expx(AGE-AGEref)fsigmoid(x,E50,λ)=xλxλ+E50λfcentral(x)=fsigmoid(x,θ12,1)fCLmaturation=fsigmoid(PMA,θ8,θ9)fQ3maturation=fsigmoid(AGE+40weeks,θ14,1)fopioids(x)=1,absence of opiatesexpx·AGE,presence of opiatesfAl-Sallami=0.88+0.121+AGE/13.4-12.79270·WGT6680+216BMI,males1.11+-0.891+AGE/7.1-1.19270·WGT8780+244BMI,femalesV1,arterial(L)=θ1fcentral(WGT)fcentral(WGTref)V1,venous(L)=V1,arterial1+θ17(1-fcentral(WGT))V2(L)=θ2WGTWGTreffageing(θ10)V3(L)=θ3fAl-SallamifAl-Sallami,reffopioids(θ13)CL(L \ min)=θ4,maleθ14,female}WGTWGTref0.75fCLmaturationfCLmaturation,reffopioids(θ11)Q2,arterial(Lmin)=θ5V2/V2,ref0.75(1+θ16(1-fQ3maturation))Q2,venous(Lmin)=Q2,arterial·θ18Q3(Lmin)=θ6V3/V3,ref0.75fQ3maturationfQ3maturation,ref

where the reference patient is male, 35 years old, 70 kg and 1.7 metres tall. All parameters values θ can be found in the original publication of [4].

Acknowledgements

We would like to thank Fredrik Bagge Carlson with JuliaHub for support in setting up the computational framework that we have used in the implementation of our proposed methodology. We would also like to thank Prof. Mats Karlsson’s research group at the Division of Pharmacometrics, Department of Pharmacy, at Uppsala University, for rewarding discussions regarding the methodology.

Author contributions

All authors designed the study, Wahlquist performed the coding, Wahlquist and Soltesz analyzed the results. All authors authored, critically reviewed and approved the manuscript for submission.

Funding

Open access funding provided by Lund University. This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. All authors are members of the ELLIIT Strategic Research Area at Lund University.

Data availibility

Code and data for reproducing the results is available in [27].

Declarations

Conflict of interest

The authors declare no conflict of interest.

1 The original data set in [4] comprises of 1033 individuals but two of them were excluded due to lack of observation data.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Sheiner LB Beal SL Evaluation of methods for estimating population pharmacokinetics parameters. i. Michaelis-Menten model: routine clinical pharmacokinetic data J Pharmacokinet Biopharm 1980 8 6 553 571 10.1007/BF01060053 7229908
Sheiner LB, Beal SL (1980) Evaluation of methods for estimating population pharmacokinetics parameters. i. Michaelis-Menten model: routine clinical pharmacokinetic data. J Pharmacokinet Biopharm 8(6):553–571. 10.1007/BF010600537229908
2. Rackauckas C, Ma Y, Noack A, et al (2020) Accelerated predictive healthcare analytics with pumas, a high performance pharmaceutical modeling and simulation platform. bioRxiv 10.1101/2020.11.28.402297
3. Marsh B White M Morton N Pharmacokinetic model driven infusion of propofol in children Brit J Anaesth 1991 67 1 41 48 10.1093/bja/67.1.41 1859758
Marsh B, White M, Morton N et al (1991) Pharmacokinetic model driven infusion of propofol in children. Brit J Anaesth 67(1):41–48. 10.1093/bja/67.1.411859758
4. Eleveld DJ Colin P Absalom AR Pharmacokinetic–pharmacodynamic model for propofol for broad application in anaesthesia and sedation Brit J Anaesth 2018 120 5 942 959 10.1016/j.bja.2018.01.018 29661412
Eleveld DJ, Colin P, Absalom AR et al (2018) Pharmacokinetic–pharmacodynamic model for propofol for broad application in anaesthesia and sedation. Brit J Anaesth 120(5):942–959. 10.1016/j.bja.2018.01.01829661412
5. Jonsson EN Karlsson MO Automated covariate model building within NONMEM Pharm Res 1998 15 9 1463 1468 10.1023/a:1011970125687 9755901
Jonsson EN, Karlsson MO (1998) Automated covariate model building within NONMEM. Pharm Res 15(9):1463–1468. 10.1023/a:10119701256879755901
6. McComb M Bies R Ramanathan M Machine learning in pharmacometrics: opportunities and challenges British J Clin Pharmacol 2022 88 4 1482 1499 10.1111/bcp.14801
McComb M, Bies R, Ramanathan M (2022) Machine learning in pharmacometrics: opportunities and challenges. British J Clin Pharmacol 88(4):1482–1499. 10.1111/bcp.14801
7. Bräm DS Parrott N Hutchinson L Introduction of an artificial neural network-based method for concentration-time predictions Pharmacometr Syst Pharmacol 2022 11 6 745 754 10.1002/psp4.12786
Bräm DS, Parrott N, Hutchinson L et al (2022) Introduction of an artificial neural network-based method for concentration-time predictions. Pharmacometr Syst Pharmacol 11(6):745–754. 10.1002/psp4.12786
8. Janssen A Hoogendoorn M Cnossen MH Application of SHAP values for inferring the optimal functional form of covariates in pharmacokinetic modeling Pharmacometr Syst Pharmacol 2022 11 8 1100 1110 10.1002/psp4.12828
Janssen A, Hoogendoorn M, Cnossen MH et al (2022) Application of SHAP values for inferring the optimal functional form of covariates in pharmacokinetic modeling. Pharmacometr Syst Pharmacol 11(8):1100–1110. 10.1002/psp4.12828
9. Sibieude E Khandelwal A Hesthaven JS Fast screening of covariates in population models empowered by machine learning J Pharmacokinet Pharmacodyn 2021 48 4 597 609 10.1007/s10928-021-09757-w 34019213
Sibieude E, Khandelwal A, Hesthaven JS et al (2021) Fast screening of covariates in population models empowered by machine learning. J Pharmacokinet Pharmacodyn 48(4):597–609. 10.1007/s10928-021-09757-w34019213
10. Sibieude E Khandelwal A Girard P Population pharmacokinetic model selection assisted by machine learning J Pharmacokin Pharmacodyn 2022 49 2 257 270 10.1007/s10928-021-09793-6
Sibieude E, Khandelwal A, Girard P et al (2022) Population pharmacokinetic model selection assisted by machine learning. J Pharmacokin Pharmacodyn 49(2):257–270. 10.1007/s10928-021-09793-6
11. Davidson JW Savic DA Walters GA Symbolic and numerical regression: experiments and applications Inf Sci 2003 150 1 95 117 10.1016/S0020-0255(02)00371-7
Davidson JW, Savic DA, Walters GA (2003) Symbolic and numerical regression: experiments and applications. Inf Sci 150(1):95–117. 10.1016/S0020-0255(02)00371-7
12. Schnider TW Minto CF Gambus PL The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers Anesthesiology 1998 88 5 1170 1182 10.1097/00000542-199805000-00006 9605675
Schnider TW, Minto CF, Gambus PL et al (1998) The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers. Anesthesiology 88(5):1170–1182. 10.1097/00000542-199805000-000069605675
13. Wahlquist Y, Bagge Carlson F, Soltesz K (2023) Fast simulation of pharmacokinetics. IFAC PapersOnline
14. Masui K Upton RN Doufas AG The performance of compartmental and physiologically based recirculatory pharmacokinetic models for propofol: a comparison using bolus, continuous, and target-controlled infusion data Anesth Analg 2010 111 2 368 379 10.1213/ANE.0b013e3181bdcf5b 19861357
Masui K, Upton RN, Doufas AG et al (2010) The performance of compartmental and physiologically based recirculatory pharmacokinetic models for propofol: a comparison using bolus, continuous, and target-controlled infusion data. Anesth Analg 111(2):368–379. 10.1213/ANE.0b013e3181bdcf5b19861357
15. Varvel JR Donoho DL Shafer SL Measuring the predictive performance of computer-controlled infusion pumps J Pharmacokinet Biopharm 1992 20 1 63 94 10.1007/BF01143186 1588504
Varvel JR, Donoho DL, Shafer SL (1992) Measuring the predictive performance of computer-controlled infusion pumps. J Pharmacokinet Biopharm 20(1):63–94. 10.1007/BF011431861588504
16. Schüttler J Kloos S Schwilden H Total intravenous anaesthesia with propofol and alfentanil by computer-assisted infusion Anaesthesia 1988 43 s1 2 7 10.1111/j.1365-2044.1988.tb09059.x 3129955
Schüttler J, Kloos S, Schwilden H et al (1988) Total intravenous anaesthesia with propofol and alfentanil by computer-assisted infusion. Anaesthesia 43(s1):2–7. 10.1111/j.1365-2044.1988.tb09059.x3129955
17. Martius G, Lampert CH (2017) Extrapolation and learning equations. In: 5th International Conference on Learning Representations, ICLR 2017, Toulon 10.48550/arXiv.1610.02995
18. Dubey SR Singh SK Chaudhuri BB Activation functions in deep learning: a comprehensive survey and benchmark Neurocomputing 2022 503 92 108 10.1016/j.neucom.2022.06.111
Dubey SR, Singh SK, Chaudhuri BB (2022) Activation functions in deep learning: a comprehensive survey and benchmark. Neurocomputing 503:92–108. 10.1016/j.neucom.2022.06.111
19. Orzechowski P, La Cava W, Moore JH (2018) Where are we now? a large benchmark study of recent symbolic regression methods. In: Proceedings of the Genetic and Evolutionary Computation Conference, Kyoto, pp 1183—1190, 10.1145/3205455.3205539
20. Sahoo SS, Lampert CH, Martius G (2018) Learning equations for extrapolation and control. In: Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholm, Sweden, Proceedings of Machine Learning Research, vol 80. PMLR, pp 4439–4447, 10.48550/arXiv.1806.07259
21. Kim S Lu PY Mukherjee S Integration of neural network-based symbolic regression in deep learning for scientific discovery IEEE Trans Neural Networks Learn Syst 2021 32 9 4166 4177 10.1109/TNNLS.2020.3017010
Kim S, Lu PY, Mukherjee S et al (2021) Integration of neural network-based symbolic regression in deep learning for scientific discovery. IEEE Trans Neural Networks Learn Syst 32(9):4166–4177. 10.1109/TNNLS.2020.3017010
22. Udrescu SM, Tan A, Feng J, et al (2020) AI feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity. In: Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS’20, pp 4860–4871, 10.48550/arXiv.2006.10782
23. Koza JR Genetic programming: on the programming of computers by means of natural selection 1992 Cambridge MIT Press
Koza JR (1992) Genetic programming: on the programming of computers by means of natural selection. MIT Press, Cambridge. 10.1007/BF00175355
24. Kingma DP, Ba J (2017) Adam: A method for stochastic optimization. In: Proceedings of the 3rd International Conference on Learning Representations (ICLR), 10.48550/arXiv.1412.6980
25. Innes M Flux: elegant machine learning with Julia J Open Source Softw 2018 10.21105/joss.00602
Innes M (2018) Flux: elegant machine learning with Julia. J Open Source Softw. 10.21105/joss.00602
26. Bezanson J Edelman A Karpinski S Julia: a fresh approach to numerical computing SIAM Rev 2017 59 1 65 98 10.1137/141000671
Bezanson J, Edelman A, Karpinski S et al (2017) Julia: a fresh approach to numerical computing. SIAM Rev 59(1):65–98. 10.1137/141000671
27. Wahlquist Y (2023) Learning pharmacometric structures. https://github.com/wahlquisty/learning-pharmacometric-covariate-structures, commit: dfafb1f
28. LeCun Y, Denker J, Solla S (1989) Optimal brain damage. In: Advances in neural information processing systems, vol 2. Morgan-Kaufmann
29. Innes M (2018) Don’t unroll adjoint: Differentiating SSA-form programs. CoRR arXiv1810.07951
30. Wahlquist Y, Morin M, Soltesz K (2022) Pharmacometric covariate modeling using symbolic regression networks. In: 2022 IEEE Conference on Control Technology and Applications (CCTA), pp 1–24, 10.1109/CCTA41146.2020.9206396
31. Mogensen PK Riseth AN Optim: a mathematical optimization package for Julia J Open Source Softw 2018 3 24 615 10.21105/joss.00615
Mogensen PK, Riseth AN (2018) Optim: a mathematical optimization package for Julia. J Open Source Softw 3(24):615. 10.21105/joss.00615
32. West GB Brown JH Enquist BJ A general model for the origin of allometric scaling laws in biology Science 1997 276 5309 122 126 10.1126/science.276.5309.122 9082983
West GB, Brown JH, Enquist BJ (1997) A general model for the origin of allometric scaling laws in biology. Science 276(5309):122–126. 10.1126/science.276.5309.1229082983
33. Anderson BJ Holford NHG Mechanism-based concepts of size and maturity in pharmacokinetics Ann Rev Pharmacol Toxicol 2008 48 303 332 10.1146/annurev.pharmtox.48.113006.094708 17914927
Anderson BJ, Holford NHG (2008) Mechanism-based concepts of size and maturity in pharmacokinetics. Ann Rev Pharmacol Toxicol 48:303–332. 10.1146/annurev.pharmtox.48.113006.09470817914927
