Learning pharmacokinetic models for in vivo glucocorticoid activation

University of Groningen

Learning pharmacokinetic models for in vivo glucocorticoid activation

Bunte, Kerstin; Smith, David J.; Chappell, Michael J.; Hassan-Smith, Zaki K.; Tomlinson,

Jeremy W.; Arlt, Wiebke; Tino, Peter

Published in:

Journal of Theoretical Biology

DOI:

10.1016/j.jtbi.2018.07.025

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from

it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date:

2018

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Bunte, K., Smith, D. J., Chappell, M. J., Hassan-Smith, Z. K., Tomlinson, J. W., Arlt, W., & Tino, P. (2018).

Learning pharmacokinetic models for in vivo glucocorticoid activation. Journal of Theoretical Biology, 455,

222-231. https://doi.org/10.1016/j.jtbi.2018.07.025

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.

ContentslistsavailableatScienceDirect

Journal

of

Theoretical

Biology

journalhomepage:www.elsevier.com/locate/jtbi

Learning

pharmacokinetic

models

for

in

vivo

glucocorticoid

activation

Kerstin

Bunte

a ,b ,∗

_,

_David

_J.

_Smith

c ,f

_,

_Michael

_J.

_Chappell

d

_,

_Zaki

_K.

_Hassan-Smith

i ,e ,g

_,

Jeremy

W.

Tomlinson

h

_,

_Wiebke

_Arlt

f ,g

_,

_Peter

_Ti

_ˇno

a ,f

a School of Computer Science, The University of Birmingham, Birmingham B15 2TT, UK

b Faculty of Science and Engineering, University of Groningen, P.O. Box 407, Groningen 9700 AK, Netherlands c School of Mathematics, The University of Birmingham, Birmingham B15 2TT, UK

d School of Engineering, University of Warwick, Coventry CV4 7AL, UK

e Departments of Endocrinology and Acute Internal Medicine, Queen Elizabeth Hospital Birmingham, Birmingham B15 2TH, UK f Institute of Metabolism and Systems Research, University of Birmingham, UK

g Centre of Endocrinology, Diabetes and Metabolism, Queen Elizabeth Hospital Birmingham, Birmingham Health Partners, UK h Oxford Centre for Diabetes, Endocrinology and Metabolism, NIHR Oxford Biomedical Research Centre, University of Oxford, Oxford, UK i Centre for Applied Biological and Exercise Science, Coventry University, Coventry, UK

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 23 October 2017 Revised 3 July 2018 Accepted 21 July 2018 Available online 23 July 2018

Keywords: Dynamic systems Pharmacokinetics Identifiability analysis Perturbation analysis 11 β-HSD activity

In vivo glucocorticoid activation Probabilistic models

Gaussian mixture model Expectation maximization Clustering

Partially observed time series analysis

a

b

s

t

r

a

c

t

Tounderstandtrendsinindividualresponsestomedication,onecantakeapurelydata-drivenmachine learningapproach,oralternativelyapplypharmacokineticscombinedwithmixed-effectsstatistical mod-elling. To take advantageof the predictivepower of machinelearning and the explanatorypower of pharmacokinetics,weproposealatentvariablemixturemodel forlearningclustersofpharmacokinetic models demonstrated on aclinical data set investigating11β-hydroxysteroid dehydrogenaseenzymes (11β-HSD)activityinhealthyadults.Theproposedstrategyautomaticallyconstructsdifferentpopulation modelsthatarenotbasedonpriorknowledgeorexperimentaldesign,butresult naturallyas mixture componentmodelsofthegloballatentvariablemixturemodel.Westudytheparameterofthe underly-ingmulti-compartmentordinarydifferentialequationmodelviaidentifiabilityanalysisontheobservable measurements,whichrevealsthemodelisstructurally locallyidentifiable.Furtherapproximationwith aperturbationtechniqueenablesefficienttrainingoftheproposedprobabilisticlatentvariable mixture clusteringtechniqueusingEstimationMaximization.Thetrainingontheclinicaldataresultsin4clusters reflectingtheprednisoneconversionrateoveraperiodof4hbasedonvenousbloodsamplestakenat 20-minintervals.Thelearnedclustersdifferinprednisoneabsorptionaswellasprednisone/prednisolone conversion.Inthediscussionsectionweincludeadetailedinvestigationoftherelationshipofthe phar-macokineticparametersofthetrainedclustermodelsforpossibleorplausiblephysiologicalexplanation andcorrelationsanalysisusingadditionalphenotypicparticipantmeasurements.

© 2018TheAuthors.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Glucocorticoids are widely used, with prescriptions for up to 2.5% of the population for immunomodulatory and anti-inflammatory effects in a number of disease states (van Staa et al., 20 0 0 ).Endogenousglucocorticoidhormones(including corti-solandcortisone)areavitalpartofnormalmetabolismand phys-iologicalfunction.Theyareproducedby theadrenalcortex under the regulation of the HPA (hypothalamic-pituitary-adrenal) axis, inadditiontoenzymatic actionintissue.Cortisolincreasesblood

∗ _{Corresponding author at: Faculty of Science and Engineering, University of}

Groningen, P.O. Box 407, Groningen 9700 AK, Netherlands.

E-mail address: kerstin.bunte@googlemail.com (K. Bunte).

sugar,functionsasanimmunesystemsuppressant,decreasesbone formationandsupportsthemetabolismoffat, proteinand carbo-hydrates. Unfortunately these hormones are also associated with adversefeatures includingcentralobesity,proximal myopathy, os-teoporosis, hypertension, insulin resistance, psychological effects andexcessiveskinchanges,whichservetoreflecttheiractionina rangeofmetabolicallyactivetissues. Theseeffectsregularlyaffect patientsreceivingexogenousglucocorticoidtreatmentbutare par-ticularly demonstratedintherareconditionofendogenous Cush-ing’s syndrome which can occur asthe result oftumours of the pituitary or adrenal gland or as a result of ectopic secretion of ACTH(Newell-Price et al., 2006 )andtheseareassociatedwith ex-cessmortality(Clayton et al., 2016; Hassan-Smith et al., 2012 ).

https://doi.org/10.1016/j.jtbi.2018.07.025

In recent years initiatives such as Horizon 2020 have been launched inorder to address thepublic health challenges ofour ageing population. In Europe, those aged >65 years made up 17 million ofthe populationin 1998,a number projectedto rise to 25 million by 2035 (Eurostat, 2008; Sakuma and Yamaguchi, 2013 ).Healthy lifeexpectancy hasnot keptup withthisincrease in longevity, with a gap between life expectancy and disability free life expectancyin theUK of9years forwomen and7years formenatthe ageof65 (Self et al., 2012 ). As aresultthere has been much research focus on the role of so-called pre-receptor metabolism of glucocorticoids via the 11

β

-hydroxysteroid dehy-drogenase enzymes(11

β

-HSD),inmetabolic conditions(including obesity and diabetes) as well as those associated with adverse ageing (osteoporosis and sarcopenia) which are similar to those foundinglucocorticoidexcess(Gathercole et al., 2013 ).

The two isozymes of 11

β

-HSD regulate glucocorticoid action at a tissue level by shuttling them between active and inac-tive forms. The type 1 enzyme (11

β

-HSD1) amplifies local tis-sue glucocorticoid levels by replacing the C11-keto group with a C11-hydroxyl group, converting endogenous cortisone to cortisol (Gathercole et al., 2013 ).Thisactivityisalsocriticallyimportantfor exogenously administered synthetic steroids, such asprednisone, whichisconvertedby11

β

-HSD1toitsactiveform(prednisolone). Thetype2enzyme11

β

-HSD2,ontheotherhand,catalyzesinvivo

mainly1_{theoppositereactionto11}

_β

_{-HSD1,enhancingthe}

inactiva-tionofcortisol/prednisolonetocortisone/prednisone.Thequestion astowhetheralteredmetabolismduetomoresubtlechanges re-latedtoageingorindividualgeneticdifferencesforexamplecould resultinchangesinglucocorticoid responsivetissuesandaccount foradverse tissue effectsis acompelling one. Assessingthe vari-ation in 11

β

-HSD1 activitybetweenindividuals,and within indi-viduals asa resultofageingandlifestylechanges,andthe result-ingmorbidityisasignificantquestioninthefieldofmetabolic re-search. 11

β

-HSD1iswidely expressed inmetabolically active tis-sues including liver, adipose, muscle, bone, skin and the central nervoussystemandtheenzymehasbeenimplicatedinthe patho-genesisofassociateddiseases(Cooper et al., 1993; Tiganescu et al., 2013 ). Cell culture and animal models have suggested that 11

β

-HSD1 is a major regulator of obesity andof the features of glu-cocorticoid excess asseen inCushing’s Syndrome (Clayton et al., 2016; Markey et al., 2016; Morgan et al., 2016a; 2016b; 2014 ). Phar-maceutical companies have developed a number of selective in-hibitorsof11

β

-HSD1andareassessingtheirtherapeuticpotential. There remains a lack of consensus on the most appropriate biomarker to measure 11

β

-HSD1 activity, which include urine steroid metabolite ratios after24 h collections,tissue biopsies to measure activity andgene expression and dynamictestssuch as theprednisolonegenerationtest.Therearelimiteddataonthe lat-ter test,whichinvolves administrationoforal prednisoneand se-rial blood testsformeasurement ofprednisone andprednisolone levels, representativeofinvivo activationofthissynthetic gluco-corticoid.Thisinformationcouldinformfuturestudyprotocols.

Prednisone hasan identicalaffinity for11

β

-HSD1ascortisone and the interconversion of oral prednisone to prednisolone has beenusedasamarkerofpredominantlyhepatic11

β

-HSD1 activ-ity (reflectingfirst pass metabolism) (Gathercole et al., 2013 ). To dateonly afew studieshaveusedtheprednisone generationtest to gain insight into the potential benefits forwell-being, healthy ageing and personalized medicine: Tomlinson et al. (2007) in-vestigated the effects of 11

β

-HSD1 inhibition in different com-partments withregard to adipose/fattissue based on serum cor-tisol and prednisolone generation in 7 healthy male volunteers;

Cooper et al. (2002) looked at in vitro (cell culture) activity (as

1 11 _{β-HSD2 is also catalyzing 11 β-HSD1 activation, but less efficient.}

opposedtoinvivoactivityaswe investigatein thiscontribution) while Hundertmark et al. (1997) and Chen et al. (1990) looked atpharmacokineticsin6healthy maleswithIVadministration of prednisoloneasopposedtoprednisone.

Invitrobiochemicalanalysisofserumprovidesamethodfor as-sessingtheactivitylevelsoftheseenzymes.Howeveritisunclear how informative serum activity data are regarding the dynamic processesoccurringinvivo. Methodsaiming toanswerthis ques-tion include in vitro-to-in vivo extrapolation (IVIVE) (Cho et al., 2014 ) techniques, whichhave become an importanttool for pre-diction of human effective dosages. However, as pointed out in

Sager et al. (2015) , IVIVE in general requires considerably more experimental and in silicodata than alternative static models. A moredirectmeasureofenzymaticactivityinvivoistointroducea prescribed doseofthe pharmacologicalcortisone analogue, pred-nisone,andtotaketimeseriesdataoftheresultingblood concen-trations,alongwiththeactivemetaboliteprednisolone.Thispaper tacklestheproblemofhowtoanalyseadatasetconsistingofsuch timeseriesfromagroupofhealthyvolunteers.

A classical pharmacokinetic approach to dealing with such datasets is to combine multi-compartment ordinary differential equation model with a mixed effects statistical model of inter-personvariation(Beal and Sheiner, 1982; Owen and Fiedler-Kelly, 2014 ). Themean(“fixedeffect”) andstandarddeviation(“random effect”) associated withthe rate constantsrepresenting reactions between prednisone and prednisolone then provide measures of central tendency and variability in 11

β

-HSD1/2 activity through thepopulationunderstudy.Forsuchsystemsthereisoftenlimited access forinputs or perturbations and the mathematical models that are generated invariably include state variables with associ-atedmodelparameterswhichareunknownandcannotbedirectly measured.Theselimitationscan causeissueswhenattempting to inferorestimate unknown modelparameters fromsets of obser-vationsandthiscan severelyhindermodelvalidation.It is there-forehighlydesirabletohaveaformalapproachtodeterminewhat additionalinputs and/or measurements are necessaryin order to reduce, orremove theselimitations andpermit the derivation of modelsthat canbe usedforpractical purposes withgreater con-fidence. Structural identifiabilityarises in the inverseproblem of inferring from the known, or assumed, properties of a systema suitablemodelstructureandestimatesforthecorrespondingrate constantsandotherparameters.Theanalysisconsidersthe unique-ness (or otherwise) ofthe unknown model parameters from the input-output structurecorresponding to proposed experimentsto collectdataforparameterestimation(underanassumptionofthe availability of perfect, noise-free observations). Thisis an impor-tant, butoftenoverlooked, theoretical prerequisite to experiment design,systemidentificationandparameterestimation,since esti-matesforunidentifiableparametersare effectivelymeaningless.If parameterestimatesaretobeusedtoinformaboutinterventionor inhibitionstrategies,orothercriticaldecisions,thenitisessential thattheparametersbeuniquelyidentifiable.Inthispapera struc-tural identifiabilityanalysisof a linearcompartmental model de-velopedtocharacteriseprednisonekineticsisperformedusingthe Laplacetransformapproach.Thisanalysisdemonstratesthat from a structural perspective the model is structurallylocally identifi-ableforthegivensystemobservations,thusprovidingmore confi-denceintheresultsobtainedforsubsequentnumericalparameter estimationusingactualtimesseriesdatafortheseobservations.

A “data driven” approachto analysing such a data set would be to apply unsupervised clustering methods from the field of machine learning such as k-means (Bishop, 1995; Lloyd, 1982 ), selforganizingmaps(Kohonen, 1982 ) andGaussianmixture mod-els (Feldman and Langberg, 2011 ), in which the data from each participant are considered as a vector from a high dimensional space. Clustering in the space of observed data either implicitly

Table 1

Subject characteristics of the 12 healthy adults recruited from the local population. The last column shows the maximum likelihood cluster membership estimate (see Section 3 ).

ID Age (years) Sex BMI (kg/m 2 ) Total Fat Mass (kg) Total Lean Mass (kg) Fasting Glucose (mmol/L) Cholesterol (mmol/L) Cls

1 20 Male 25.0 16.1 59.0 5.1 4.2 C2 2 69 Male 24.0 15.3 45.5 4.5 5.1 C4 3 26 Male 25.9 17.5 57.7 5.1 4.3 C2 4 57 Male 27.5 18.5 61.7 5.1 5.8 C1 5 25 Male 24.1 20.2 46.1 5.5 3.4 C1 6 54 Male 25.0 16.5 45.0 4.8 5.1 C2 7 23 Female 22.2 17.2 40.7 4.3 4.3 C3 8 64 Female 22.0 18.9 39.8 4.9 4.4 C4 9 24 Female 21.4 18.2 42.6 4.5 5.1 C4 10 50 Female 24.7 23.3 35.8 4.1 5.7 C2 11 20 Female 19.9 10.9 36.3 5.0 3.8 C4 12 60 Female 29.1 33.0 44.8 4.5 6.5 C2

or explicitly assumes a certain metric or structure of the data. In this contribution we will pursue a hybrid multidisciplinary approach of model based clustering integrated with structural identifiabilityanalysisforinterpretationofparameterrelationsand perturbationtheory to reduce the dimensionality ofthe parame-ter space.The core of ourcontribution is based on interpretable probabilisticinferential models,aiming atgroupingindividuals in the space of pharmacokinetic models based on observed data. Eachgroup/clusteristhereforerepresentedbyaprototypical prob-abilisticmodel with a specific pharmacokineticparameterisation. This way our proposed strategy automatically constructs differ-ent“population” models,that are thereforenot definedbasedon prior knowledge or experimental design, butcome out naturally asmixture component models of the globallatent variable mix-ture model. In contrast to data driven clustering techniques, we can analyse the parameter relationships and investigate possible orplausiblephysiologicalexplanation.The investigationoffurther phenotypicmeasurementsofindividualsmoreprobabletobe rep-resentedbythesameclustermodelmightleadtonewhypothesis ofinterestingbiomarkersforfutureinvestigationandclinical stud-ies. This contribution therefore reveals both the capabilities and limitations ofpharmacokinetic modelling combined with param-eterestimationandmachine learning, demonstratedon aclinical datasetforprednisone conversionasan exampleofits potential broaderapplicationtomodellingofinvivobiochemicalsystemsin heterogeneouspopulations.

2. Materialsandmethods

2.1.Clinicaldata

Theinvestigationswereperformedin12healthyadults(6men and6 women)recruited fromthe local population atthe Queen Elizabeth Hospital Birmingham with subject characteristics sum-marized in Table 1 . Inclusion criteria included body mass index between20 to 30kg/m2_{, femalesin the follicular phase of their}

menstrual cycle, and post-menopausal subjects off estrogen re-placementtherapy. Exclusion criteria included pregnancy, signifi-cantpast medicalhistory(like diabetesmellitus),ischaemicheart disease,cerebrovascular disease,respiratory disease andepilepsy, useofdrugsincludingglucocorticoids,beta-blockers,dopamine ag-onistsandanticoagulants.AclinicalstudyinBirminghamwas car-riedout betweenOctober 2010 and March 2013. Participants ar-rived at the NIHR-Wellcome Trust Clinical Research Facility in a fasted state by 8:30 AM. Baseline blood tests were takenat ap-proximately9:00AMandanalysedforureaandelectrolytes,lipids, glucose(RocheModularSystem),insulin(colourimatricELISAfrom Mercodia)in additionTSH andfree T4(Advia Centaur; Bayer Di-agnostics) were sent. 10 mg of prednisone was then adminis-tered orally with additional venous blood samples taken at

20-Fig. 1. Three-component model schematic. Fast processes ( P  L ) are represented with bold arrows.

min intervals over a period of 4 h with serum extracted and analysed for cortisol and cortisone serum concentrations by liq-uid chromatography-mass spectrometry as previously described (Hassan-Smith et al., 2015 ). Observationsincludingheight,weight andbloodpressure wererecordedandbody composition was as-sessed using Dual-energy X-ray absorptiometry (DXA) scannning (HologicDiscovery:versionApex3.0,HologicInc).

Ethical Approval. The study was approved by the Coventry and Warwickshire Research Ethics Committee (REC reference no. 07/H1211/68)andtheScientific CommitteeoftheNIHR-Wellcome Trust Clinical Research Facility at the Queen Elizabeth Hospital Birmingham.

2.2. Linearkinetics,threecompartmentmodel

The model (Fig. 1 ) consistsof three compartments, an unob-served stomach compartment S in which the prednisone formu-lation isinitially deposited afteroral ingestion, acompartment P

representingbloodconcentration(nmol/L)oftheinactive metabo-liteprednisone,andacompartment Lrepresentingblood concen-tration (nmol/L) of the active metabolite prednisolone. Reactions betweenthe three compartmentsare assumed to havelinear ki-netics.Thiscontributionfocusesonthein-depthanalysisofa lin-earkineticexampleasaproofofconcept,pavingthewayfor anal-ysiswithincreasedcomplexity,suchasnon-linear models,inthe future.Rate constantsforthe linearkinetics include:k_abs for ab-sorptionoforalprednisoneintotheblood,kPLandkLPrepresenting

11

β

-HSD1and11

β

-HSD2activityconvertingtheinactive metabo-litetoactiveform,andviceversa,andexcretionconstantskPexand kLex quantifyingexcretionfromthebloodofprednisoneand

pred-nisolonerespectively.Theenzymeactivityisassumedtobe signif-icantly fasterthaneitherthe absorptionorexcretion respectively, justifiedby thealmost immediatepresence intheblood of pred-nisoloneobservedexperimentally.

Withall reactionsassumedtohavelinearkinetics,the mathe-maticalmodelthereforetakestheform:

dS dt =−kabsS, dP dt =kabsS−

(

kPex+kPL

)

P+kLPL, dL dt =kPLP−

(

kLex+kLP

)

L,

⎫

⎪

⎪

⎪

⎪

⎬

⎪

⎪

⎪

⎪

⎭

(1)

withinitial condition[S,P,L]=[S0,0,0] attimet=0with

ob-servations:



_y 1 y2 y3



=



_{0 0 0} 0 1 0 0 0 1





_S P L



. (2)

Thismodelhassixparameters,S0,kabs,kPL,kLP,kPex andkLex and

canbewrittenmorecompactlyintheform:

d

μ

dt =A

μ

,

μ

(

0

)

=[S0,0,0]

T_{, (3)}

y=C

μ

(4)

where

μ

(

t

)

=[S

(

t

)

,P

(

t

)

,L

(

t

)

]T_{, yincludes theobservable} dimen-sionsy2 andy3(correspondingtoPandL)asselectedbymatrixC

asgiveninEq. (2) andthematrixAisdefinedby:

A=



_−k abs 0 0 kabs −

(

kPex+kPL

)

kLP 0 kPL −

(

kLex+kLP

)



. (5)

The formulation Eqs. (2) –(5) explicitly defines model parameter and model output structure, which will be analysed in detail in thefollowingsection.

2.3. Structuralidentifiabilityanalysis

Inordertoestimatethe(unknown)modelparametersfromthe dataavailableitisnecessarytoincludeinthemodeloutput struc-ture,whichcorrespondstothefunctionofthemodelvariablesthat istobe comparedwiththedata.Beforeactuallycollecting exper-imentaldataitisnecessarytotestthosemodelvariableswith re-spect to thisoutput structure foruniqueness, since estimates for unidentifiableparametersaremeaningless.Such astructural iden-tifiability analysis (Bellman and ˚Aström, 1970 ) assesses whether theobservedmodeloutputcontainsenoughinformationto deter-mine all of the model parameters uniquely (Jacquez, 1996 ), and relates onlyto the structure of themodel andoutput. Forlinear systems there are manywell-establishedtechniques for perform-inga structuralidentifiabilityanalysis(forfurtherdetails,and de-tailsofnonlinearapproaches,seethetutorialbyGodfrey and DiS- tefano III (1987) andotherworksinthesamevolumeandthebook byWalter, 1982 ).

Here the uniqueness, of the unknown parameters in a gen-eralsystemsmodelisconsidered withrespecttotheoutputs. Let

p_∈



_{⊂ R}r _{denote a vector comprisingthe unknown parameters} in the model, which belongs to an open set of admissible vec-tors(Evans et al., 2002 ).Tomaketheparameterdependenceofthe modeloutputsmoreexplicititiswritteny(t,p).

Two parameter vectors p, p∈



are indistinguishable, written

p∼ p,iftheygiverisetoidenticaloutputs:

y

(

t,p

)

=y

(

t,p

)

forallt≥ 0 .

Forgenericp∈



,theparameterpiislocallyidentifiableifthereis aneighbourhood,N,ofpsuchthat

p∈N, p∼ p impliesthat pi=pi .

In particular, ifN=



in the above definitionthen pi is globally identifiable,otherwiseitisnon-uniquely(locally)identifiable.Notice

that,for agiven output,a locally identifiableparameter cantake anyofadistinct(countable)setofvalues.Iftheredoesnotexista suitableneighbourhoodN thenpiisunidentifiableand,foragiven output,cantakean(uncountably)infinitesetofvalues.

Asystemmodelisstructurallygloballyidentifiable(SGI)ifall pa-rametersare globally identifiable;it is structurally locally identifi-able(SLI)ifallparametersarelocallyidentifiableandatleastone isnon-uniquelyidentifiable;andthemodelisstructurally unidenti-fiable(SU)ifatleastoneparameterisunidentifiable.

Anestablishedapproachtoidentifiabilityanalysisoflinear sys-tems is to take the Laplace transform, reducing the initial value problemto analgebraicone.Denoting theLaplacetransformofμ as, ¯

μ

(

s

)

=  _∞ 0 e−st

μ

(

t

)

dt, (6)

theinitialvalueproblem(3) istransformedto,

−

μ

(

0

)

+s

μ

¯

(

s

)

=A

μ

¯

(

s

)

, (7)

hence,

¯

μ

(

s

)

=

(

sI− A

)

−1

_μ

₍

₀

₎

_{. (8)}

Definingthecharacteristicpolynomial,

χ

(

s

)

=det

(

sI− A

)

=

(

s+kabs

)

·

(

s2₊

₍

_k

Pex+kPL+kLex+kLP

)

s+kPexkLP+kLexkLP+kPexkLex

)

(9)

thesolutioninLaplacespacefortheobservablecomponentsis,

ˆ P

(

s

)

ˆ L

(

s

)

=

ˆ

μ

2

(

s

)

ˆ

μ

3

(

s

)

=S0kabs

χ

(

s

)

kLP+kLex+s kPL

. (10)

Expressingthesolutionintheformofrationalfunctionsyields:

ˆ

μ

2

(

s

)

=



1 s+



2 s3₊



3s2+



4 s+



5 (11) ˆ

μ

3

(

s

)

=



6 s3₊



3s2+



4 s+



5 , (12)

wherethe coefficientsof the powers ofs in thenumerators and denominators ofEqs. (11) and(12) are termed the “moment in-variants” forthese input/outputexpressions and, interms of the originalmodelparameters,aregivenby:



1=S0kabs, (13)



2=S0kabs

(

kLP+kLex

)

, (14)



3=kabs+kPex+kPL+kLex+kLP, (15)



4 =kabs

(

kPex+kPL+kLex+kLP

)

+kPexkLP+kLexkLP+kPexkLex, (16)



5=kabs

(

kPexkLP+kLexkLP+kPexkLex

)

, (17)



6=S0kabskPL. (18)

Themomentinvariantsforthesystemareassumedtobe measur-able(known)throughtheobservations,andareconsideredunique. Thesystemistermed(globally/locally)identifiableifthemapping



:[S0,kabs,kPL,kLP,kPex,kLex] →[



1,



2,



3,



4,



5,



6] (19)

is(globally/locally)invertible.OursystemEq. (1) islocally identifi-ablewiththreepossiblesolutionsasfollowsfromthesixmoment invariantsEqs. (13) –(18) :Eq. (13) impliesthattheproduct

isstructurallygloballyidentifiable(SGI).WithEq. (20) substitutedin

Eq. (14) itfollowsthat



2=S0kabs

(

kLP+kLex

)

and:



2=



1

(

kLP+kLex

)

⇒ thesumkLP+kLex=



_

2 1

(21)

isSGI.ThissubstitutedintoEq. (15)



3=kabs+kPex+kPL+

(

kLP+kLex

)

leadsto:



3=kabs+kPex+kPL+



_

2 1

implyingthesum

kabs+kPex+kPL=



3−



_

2 1

(22)

isSGI.Furthermore,substitutingEq. (20) inEq. (18) yields:



6=S0kabskPL=



1kPL

⇒kPL=



_

6 1

(23)

isSGI.NotealsofromEqs. (23) and(22) thatthesum

kabs+kPex=



3−



_

2 1−



6



1 (24) isSGI.Eq. (17)



5 kabs

=kPexkLP+kLexkLP+kPexkLex

substitutedintoEq. (16) leadsto:



4=kabs

(

kPex+kPL+kLex+kLP

)

+_k



5

abs

andsubstitutingEqs. (21) and(22) resultsin:



4=kabs



−kabs+



3−



_

2 1 +



2



1



+



5 kabs ,

whichyieldsacubicequationink_abs:

⇒k3

abs−



3k2abs+



4kabs−



5=0 . (25)

FromanalysisusingDescartes’ruleofsignsitcanreadilybeshown thatEq. (25) hasthree changesofsignsinthecoefficients, which meansithasmaximalthreepositive(real)roots.Sincethenegative polynomial f

(−k

abs

)

hasnochangeofsignithasnonegativeroots.

Furthermore,Eq. (17) canbewrittenas



5=kabs

(

kPex

(

kLP+kLex

)

+kLexkLP

)

. (26)

Insummarythestructuralidentifiabilityanalysisyields:

(1) kabs isstructurallylocallyidentifiable(SLI)withupto3possible

solutions.

(2) S0isSLI(basedonEq. (20) ).

(3) kPex isSLI(followsfromEq. (24) ).

(4) kLex andkLPareSLI(sincekLex+kLPisSLIseenbyEq. (21) and kLexkLPisSLIbecauseof(1),(3),Eqs. (21) and(26) ).

(5) kPLisSGI(seeEq. (23) ). 2.4.Dimensionalanalysis

Because the model is linear, the dependent variables may be scaled arbitrarily; choosing the initial stomach concentration S0

as the scaling factor and denoting dimensionless variables with primeswehave,

S=S0S, P=S0P, L=S0L. (27)

Takingastime-scaletheinverseoftherateofabsorptionfromthe stomach,i.e.t=k−1_abst,theproblemcanbewrittenintermsof di-mensionlessvariablesas,

dS dt =−S, dP dt =S−

η

P+

ρ





P+1



L, dL dt =

ρ



P−

η

L+ 1





L,

⎫

⎪

⎪

⎪

⎪

⎬

⎪

⎪

⎪

⎪

⎭

(28)

with initial condition

(

S,P,L

)

=

(

1,0,0

)

 at time t=0. The groups

η

P=kPex/kabs and

η

L=kLex/kabs arethedimensionless

ex-cretion rates of prednisone and prednisolone respectively,

ρ

=

kPL/kLP is the ratio between the rates of forward and backward

conversion betweenthe inactive and active metabolites, and



=

k_abs/kLPtheratiobetweentherateofstomachabsorptionandrate

ofbackwardconversionfromtheactivetoinactivemetabolite.The assumptionthatconversionbetweenmetabolitesoccursonafaster time-scalethanabsorptionandexcretionimpliesthat



_1.

2.5. Quasi-steadyapproximation

Whilemodel(28) islinearandthereforecanbe readilysolved viamatrixexponentials, itispossibletoexploit thesmall param-eter



to yield an even simplersystem with two fewer free pa-rameters. Eq. (28) for S decouples and has the analytic solution

S

(

t

)

=e−t. Substituting for S and adding the remaining equa-tions yields thefollowing equation forthe dynamicsof the total inactiveandactivemetabolitesintheblood:

d

(

P+L

)

dt =e−t



−

η

PP−

η

LL, (29)

withinitialconditionP+L=0att=0.

Eq. (29) isexact. The presenceof the smallparameter



moti-vates seeking an approximatesolution witha smaller numberof parameters. Examining Eq. (28) and retaining only terms O(1/



) yieldsthequasi-steadyapproximation,

ρ

P≈ L, (30)

therefore, P+L≈

(

1+

ρ

)

P≈

(

1+

ρ

)

ρ

−1_{L. Intuitively, this} rela-tion can be interpreted as 11

β

-HSD1 activity being sufficiently rapidthattheratiobetweenactiveandinactivemetabolitesis ap-proximatelyconstantoverthetime-scalesassociatedwith absorp-tionandexcretion.SubstitutingintoEq. (29) wethenhavethe ap-proximatemodelfortotalmetabolitedynamics,

d

(

P+L

)

dt ≈ e −t −

ηP

+

ρη

L 1+

ρ



(

P+L

)

, (31)

withinitial conditionP₊L₌0att₌0.Theanalytic solutionis,

(

P+L

)(

t

)

≈

_η

1+

ρ

P− 1+

ρ

(

η

L− 1

)

e−t− exp

−

η

P+

ρη

L 1+

ρ

t



. (32)

Thedimensionlessmetaboliteconcentrationscanthenbe deter-minedbysubstitutingexpression(30) into(32) togive

P

(

t

)

≈

_η

1 P− 1+

ρ

(

η

L− 1

)

e−t− exp

−

η

P+

ρη

L 1+

ρ

t



, (33) L

(

t

)

≈

ρ

η

P− 1+

ρ

(

η

L− 1

)

e−t− exp

−

η

P+

ρη

L 1+

ρ

t



. (34)

The above can be derived formally as the leading order terms ina perturbation expansion P

(

t

)

=P₀

(

t

)

+



P₁

(

t

)

+..., L

(

t

)

=

Fig. 2. Comparison between numerical solution of the full six-parameter model

(1) and the approximate solution with four parameters given by Eqs. (36) and (37) , based on parameter values k abs = 4 , k PL = 90 , k LP = 30 , k Pex = 2 , k Lex = 1 . 6 , S 0 = 10 . Following re-dimensionalisation, the approximate solution of thesystemEq. (1) can beexpressedintermsoftheoriginal vari-ablesandparametersas,

S

(

t

)

≈ S0exp

(

−kabst

)

, (35)

P

(

t

)

≈ S0kabskLP

kPLkLex+kLPkPex− kabs

(

kPL+kLP

)

·



exp

(

−kabst

)

− exp



−kLPkPex+kPLkLex kPL+kLP t



, (36) L

(

t

)

≈ S0kabskPL

kPLkLex+kLPkPex− kabs

(

kPL+kLP

)

·



exp

(

−kabst

)

− exp



−kLPkPex+kPLkLex kPL+kLP t



. (37)

The following four-parameter approximate model μ(t; w) can thenbedefinedfortheobservedvariablesP≈ μ2,L≈ μ3:

μ

2

(

t; w

)

=w2

(

exp

(

−w1t

)

− exp

(

−w1w3t

)

)

, (38)

μ

3

(

t; w

)

=w2w4

(

exp

(

−w1t

)

− exp

(

−w1w3t

)

)

, (39)

wheretheparametervectorw₌[w1,w2,w3,w4]isrelatedtothe

physicalparametersby,

w1=kabs (40)

w2=

S0kabskLP

kPLkLex+kLPkPex− kabs

(

kPL+kLP

)

(41) w3= kLPkPex+kPLkLex kabs

(

kPL+kLP

)

(42) w4= kPL kLP. (43)

InFig. 2 we seeexcellent agreementbetweenthenumerical solu-tion ofthe six-parametermodelEq. (1) computedwiththe Mat-lab function ode15s and the four parameters approximate solu-tion givenby Eqs. (36) and(37) ,withparameter values k_abs=4, kPL=90,kLP=30,kPex=2,kLex=1.6,S0=10yielding



=0.133.

2.6.Maximumlikelihoodestimationforclusteringof pharmokokineticmodels

The data for the absorption of prednisone and conversion to prednisolone (abbreviated by P and L respectively) exhibit large variationsinhealthyindividuals.Therefore,thequestionarisesas to whether there are groups of people with similar trajectories over the course of time and if those groups havecommon phe-notypicalcharacteristics.

We assume the generating process for the observed data can be reasonably well approximated by the simplified model

Eqs. (38) and (39) parameterised by w. The aim is to find pa-rametervectorswkcorrespondingtodifferentmodelsexplainingk groupsofconversionbehaviouroverthecourse oftime,suchthat prednisone (P) andprednisolone (L) measurements in the serum ofsubjectswithinaclusteraremoresimilarcomparedtoanother cluster. We use Maximum Likelihood Estimation of a Gaussian MixtureModelexplainedinthefollowing.ThedatasetD=

{

Dj

_}

N

j=1

iscomposedof Nindividuals.Foreach participantj we assume a collection Dj₌

_{

_yj

c

(

tmj

)

}

o_cj

m=1 of measurements of P and L attime

pointst_mj. Here y_cj with c∈{2, 3} corresponds to noisy measure-mentsofPandL,whicharemodelledbycomponentsμcofthe dy-namicalsystemMapproximatedbyEqs. (38) and(39) .Thenumber ofavailableobservationsofsubstancecandsubjectjisdenotedby

ocj andmayvaryacrosssubjects.Weassumethattheobservations fromsubstancescarenormallydistributedwithmeancomponent μcandsubstancerelatedvariance

σ

c2.Therefore,theprobabilityof a measurement ycj

(

tmj

)

being produced by component

μ

c

(

tmj,wk

)

is

P

(

ycj

(

tmj

)

|

wk

)

=N

(

ycj

(

tmj

)

;

μ

c

(

tmj,wk

)

,

σ

c2

)

(44) Fig. 3 illustrates our machine learning model and assumptions on the data. The curves denote the cluster model concentra-tions of prednisone μ2(t; wk) (red) and prednisolone μ3(t; wk) (blue), whereas differently shaped markers represent individual measurements from participants for both compounds at differ-ent time points y_cj

(

t_mj

)

respectively. The solid black line illus-trates the normally distributed observations y_cj₌₃

(

200

)

centered at

μ

c=3

(

200,wk

)

withvariance

σ

c2=3. Since measurements of

dif-ferent individualsare considered independent we write the like-lihood of the data given parameters

₌

{{

w_k,P

(

k

)

}

K

k=1,

{

σ

c

}}

as P

(

D

|

)

=N

j=1P˜

(

Dj

|

)

.Wemaximizetheloglikelihood:

logP

(

D

|

)

= N 

j=1

logP˜

(

Dj

|

₎

₍₄₅₎

usingamixturemodel

˜ P

(

Dj

|

)

= K  k=1 Pk· ˜P

(

Dj

|

k

)

(46) withpriorsPkand

˜ P

(

Dj

|

k

)

= 3  c=2



_oj c  m=1 P

(

ycj

(

tmj

)

|

wk

)



1/oj c . (47)

Weweighteachobservationtoaccountforthepotentiallyvarying numberofsuccessfulmeasurementsocjperparticipant.

FortheEMalgorithmwe assumethatthelatentBernoulli ran-dom indicator variable Z_kj is 1 ifthe data collectionDj _was gen-eratedfromthemodelparameterisedbywk and0otherwise.The completelog-likelihoodisgivenby

L

(

{

Z_kj

}

,D

)

= N  j=1 K  k=1 Z_kjlog

(

P˜

(

Dj

_|

k

)

Pk

)

. (48)

Fig. 3. Illustration of the data, variables and assumptions to train our cluster algorithm. P and L represent the compartments for blood concentration of prednisone and prednisolone respectively. The markers show example data from subjects with ID 1, 3 and 10 (see Table 1 ).

The maximum likelihood estimate (MLE) is determined by marginalizingthelikelihoodoftheobserveddatabyiteratingover 2steps:

E:Sincethecompletelikelihoodisnotknownwecalculatethe expectedvalueoftheloglikelihoodfunctionwithrespectto the conditional distribution of Z givenD giventhe current estimateoftheparameters

(t)_:

Q

(

|

(t)

₎

₌N j=1 K  k=1 E[Z_kj]· log

(

P˜

(

Dj

_|

k

)

Pk

)

Q

(

|

(t)

)

= N  j=1 K  k=1

γ

(t) jk · log

(

P˜

(

D j

_|

k

)

Pk

)

with (49)

γ

(t) jk = ˜ P

(

Dj

|

(t) k

)

P( t) k  lP˜

(

Dj

|

( t) l

)

P (t) l . (50)

M:Findtheparametersmaximizingthefollowingquantity:

(t+1) _=argmax

Q

(

|

(t)

₎

₍₅₁₎

3. Results

We performed experiments varying the number of clusters from2to6possiblemodelstorepresentthetimecourseof pred-nisone and prednisolone for the respective number of groups of subjects.Foreach numberofclusterswe useleave-one-out cross validationofthe12peoplewith5independentrepetitions result-ingin60clusteringsforeachexperiment. Thealgorithmisalways initializedwithindividualfitsofrandomlychosentrainingsubjects corresponding to the number of clusters assumed. The resulting medianlog-likelihoodEq. (48) and50%IQRversusthenumberof clustersisshownintheboxplotinFig. 4 .Incrementingthenumber ofclusters up to 4 increases the performance considerably, after thataddinganotherclusterdoesnotimprovetheclustering signif-icantly.Basedonthis“elbowcriterion” (Ketchen and Shook, 1996 ) wefurtheranalyze4classclusteringresultsindetail.

For the investigation of the models and parameters we first computetheprobabilitiesPijbycountinghowoftentwosubjectsi andjappearedtogetherinthesameclusterthroughoutthe60 in-dependentruns,assuming4clusters.Wedefine1_{− Pi j}asthe pair-wise distanceand perform completelinkage clusteringas shown inthedendrogramFig. 5 .Thecut-off valueof4clustersisusedas clusterassignmentZ_kjtoinitializeourEMprocedureoncemoreto findthefinal model parameters fordetailedinvestigation.In this last experiment we initialize by a least squares individual group

Fig. 4. Median log-likelihood and 50% IQR boxplot versus the number of clusters used for the “elbow” method to determine the number of clusters to model the data.

Fig. 5. Complete linkage clustering based on pairwise probabilities P ij of partici-

pant ID i and ID j appearing together in the same cluster throughout the experiment assuming 4 models.

fitbasedon5repetitionsofrandomparametersandtheresulting finalmodelsareshowninFig. 6 .

4. Discussion

Thefinaltrainedclustermodelsareveryrobustwithrespectto therandominitializationinthetrainingbasedonthegiveninitial 4 cluster assignment of participants extracted from the dendro-gram Fig. 5 .Therefore,only the averageparameters cluster mod-elsareshowninFig. 6 .Cluster1resemblestheslowestabsorption of prednisone and also the least growing concentration of pred-nisoloneinthe blood.The secondslowestabsorption rateis cap-turedbycluster2.Peopleassignedtothatclusterreachhigher lev-elsofprednisolone.Peopleinthefourthclusterreachhigher con-centrations of prednisolone 20 min afteradministration of pred-nisone andafter 100 min it startsto decrease. The third cluster containingonlysubjectID7resemblesanoutlierfromthedataset.

Table 2

Relationship of the original parameters for the 4 cluster models shown in Fig. 6 . We show means and standard deviations based on 5 random initializations.

kLex = a · k Pex + b C S0 · k abs = kPL = a b 1 3.74362 (0.00976) 5.02783 (0.00614) ·k LP −0 . 19889 (0.0 0 024) 0.00480 (0.0 0 0 05) 2 5.47217 (0.00744) 6.00644 (0.00293) ·k LP −0 . 16649 (0.0 0 0 08) 0.00391 (0.0 0 0 05) 4 21.34430 (0.01347) 7.40155 (0.00212) ·k LP −0 . 13511 (0.0 0 0 04) 0.01675 (0.00999) 3 34.35457 (0.02879) 12.56240 (0.00352) ·k LP −0 . 07960 (0.0 0 0 02) 0.01116 (0.0 0 021)

Fig. 6. Four learned models of prednisone (P; dashed lines) and prednisolone (L; solid lines) blood concentrations and the actual measurements of the 12 subjects ( Table 1 ).

The conversionpatternof thissubjectisvery differentcompared toallotherindividualsinthisdataset.

With the knowledge aboutthe functional relationshipsof the parametersfromtheidentifiabilityanalysiswecaninvestigatethe trained parameters and investigate possible explanations for the behaviour.Inordertoviewtherelationshipfortheoriginal param-eters(Fig. 1 )wesolvethelinearequationsystemEqs. (40) –(43) for each ofthe cluster models. The mean andstandard deviationof thoserelationshipsforeachclusterbasedonthe5random initial-izationsisshowninTable 2 .Verysmallstandardvariations quan-tifytheabovestatementthatthefinalclustermodelsarevery ro-bustwithrespecttorandominitialization.Theexcretionsof pred-nisolonekLexandprednisonekPex exhibitalinearantiproportional

relationship: themore prednisone is excreted thelower the rate for prednisolone. Only a small interval of kPex values is possibly

dependent onthevalueofb wherekLex islarger(see Fig. 7 ).The

clustersC1,C2,C4toC3exhibitincreasingabsorptionrateS0· kabs

of prednisoneinto the blood ascan be seen fromTable 2 . Pred-nisone is converted toprednisolone by a factor of5 timesfaster thanthatforcluster1.Fortheotherclustersthisfactoris increas-ing.

Afterestimatingthefinalclusterassignmentwegotbacktoour medical expert asking for additional information about the par-ticipantsto investigatecharacteristics within thesame cluster. Of course,duetothesmallsamplesizethefollowinganalysisisonly exploratoryandwe cannotmakestrongclaims, butthe investiga-tionmightleadtonewhypothesisandprovidesomeevidencefor interestingmeasurements forfutureclinicalstudies.Therefore,we compute theAreaUndertheCurve(AUC)ofprednisolone

conver-Fig. 7. Linear relationship of k Pex and k Lex for clusters C1 − C4 from Fig. 6 . The cir-

cle mark the point were k Pex = k Lex .

sionofeachclustermodelk:

AUC=  t=240

t=0

μ

3

(

t; wk

)

(52)

andweight it with the trained mixture component

γ

jk foreach participant j resulting in an individual AUC value dependent on themembershipin eachcluster. We testthe associationbetween theseAUCvaluesandmeasurementoftheadditionalclinical satel-lite data for each individual. First we observe that there are no strongcorrelationsbetweenthe ageorBMI ingeneral, witha p -value of .45 and .128 respectively. Cortisol on the other hand is correlatedtotheAUC andexhibitsap-value of.006,whichisnot verysuprising sinceprednisolone isa syntheticderivative of cor-tisolprocessedby the sameenzymesasmodeled here. Since the

hormonalactivity of female andmale population isdifferent we furthermoreinvestigatethecorrelationdividingthegroupsby gen-der.Thereissome indicationthatspecificmeasurementscouldbe interestingassubjectforfurtherinvestigation,forexamplewe ob-serveap-valueof.029fortheBMIwithin men,whileitis.47for thefemaleparticipants.The investigationofthe phenotypic mea-surementsofindividualsthataremoreprobabletoberepresented bythesameclustermodelmightthereforeleadtonewhypothesis of potentially interesting biomarkers for future investigation and clinicalstudies.

5. Conclusionsandfuturework

The core of our contribution is based on interpretable prob-abilistic inferential models aiming at clustering pharmacokinetic modelsfortheabsorptionofprednisoneintheblood. The collec-tionofallindividualbloodmeasurementsismodelledasa proba-bilisticlatentvariablemodelwithpharmacokineticmodelsplaying therole ofmixture components.Therefore,each group orcluster isrepresentedbyaprototypicalprobabilisticmodelwithaspecific pharmacokinetic parameterisation. This way our proposed strat-egy automatically constructs different “population” models, that arethereforenotdefinedbasedonpriorknowledgeor experimen-taldesign,butcome out naturallyasmixturecomponentmodels ofthe globallatent variable mixture model.In contrast to solely data-driven clustering techniques, we can analyse the parameter relationshipsandinvestigatepossibleorplausiblephysiological ex-planation.Thestrategyissuitableforsparsemeasurements,which isespecially beneficialifthesearecollectedbyan invasive proce-dure.Ourapproach isdesignedfortimeseriesmeasurements po-tentiallytakenatdifferenttime points andisdemonstrated ona clinicaldataset investigatingtheinvivoglucocorticoid activation by11

β

-HSD1/2activityinhealthyadults.

The model was thoroughly studied by identifiability analysis andthenapproximatedusingtheperturbationmethod.Thelatent variablemixture ofpharmacokinetic models istrained by an Ex-pectationMaximization strategy, which isa widely used efficient naturalchoice for the estimation of such latent variablemodels. We achieved robust results for 4 prototypical cluster models re-semblingtheprednisone/prednisolone concentration in theblood overthecourseof240minafteradmissionofthedrugforthe12 subjects in the data set. We observed a weak correlation of the AUC of prednisolone concentration in the blood with respect to theclustermodelsandtheBMIofmalesuspects,whichdoesnot seemtobeimmanentforthefemaleparticipants.Theinvestigation offurtherphenotypicmeasurements ofindividualsmoreprobable tobe represented by the same cluster modelmight lead to new hypothesis of interesting biomarkers for future investigation and clinical studies. With the availability of more data in the future theapproachcanbeextendedtonon-linearpharmacokinetic mod-els,whilethiscontributionserves asa proofofconcept.Our pro-posalemphasise the potential of exploratoryanalysis ofpartially observedtime seriesdata by combiningpharmacokinetics, struc-turalidentifiabilityanalysis,dimensionalanalysis/perturbation the-orywithmachinelearning.

Acknowledgements

This work has received funding from the European Union’s Horizon H2020-EU.1.3.2.-Nurturing excellence by means ofcross-borderandcross-sector mobility, MSCA-IF-2014-EF-MarieSkłodowska-Curie Indi-vidual Fellowship (IF-EF) undergrantagreement no. 659104. And the authors acknowledge sup-portfromEngineeringandPhysicalSciences ResearchCouncil EP-SRCgrantnumbersEP/L000296/1andEP/N021096/1.

Furthermore we thank Theresa Brady, Pamela Jones, Claire Brown(NationalInstituteforHealthResearch-WellcomeTrust Clin-ical Research Facility, Birmingham).The clinical visit was carried outattheNationalInstituteforHealthResearch(NIHR)/Wellcome Trust Birmingham Clinical ResearchFacility. The views expressed are thoseofthe authors(s)andnot necessarilythose oftheNHS, theNIHRortheDepartmentofHealth.

References

Beal, S.L. , Sheiner, L.B. , 1982. Estimating population pharmacokinetics. CRC Critical Rev. Biomed. Eng. 8, 195–222 .

Bellman, R., ˚Aström, K., 1970. On structural identifiability. Math Biosci 7 (3), 329– 339. doi: 10.1016/0025-5564(70)90132-X .

Bishop, C.M. , 1995. Neural Networks for Pattern Recognition. Oxford University Press, Inc., New York, NY, USA .

Chen, M.F. , Shimada, F. , Kato, H. , Yano, S. , Kanaoka, M. , 1990. Effect of glycyrrhizin on the pharmacokinetics of prednisolone following low dosage of prednisolone hemisuccinate. Endocrinol. Jpn. 37 (337), 331–341 .

Cho, H.J., Kim, J.E., Kim, D.D., Yoon, I.S., 2014. In vitro-in vivo extrapolation (IVIVE) for predicting human intestinal absorption and first-pass elimination of drugs: principles and applications. Drug Dev. Ind. Pharm. 40 (8), 989–998. doi: 10.3109/ 03639045.2013.831439 .

Clayton, R., Jones, P., Reulen, R., Stewart, P., Hassan-Smith, Z., Ntali, G., Karavitaki, N., Dekkers, O., Pereira, A., Bolland, M., Holdaway, I., Lindholm, J., 2016. Mortality in patients with Cushing’s disease more than 10 years after remission: a multi- centre, multinational, retrospective cohort study. Lancet Diabetes Endocrinol. 4 (7), 569–576. doi: 10.1016/S2213-8587(16)30 0 05-5 .

Cooper, C., O’Neill, T., Silman, A., 1993. The epidemiology of vertebral fractures. European vertebral osteoporosis study group. Bone 14, 89–97. doi: 10.1016/ 8756- 3282(93)90358- H .

Cooper, M.S., Rabbitt, E.H., Goddard, P.E., Bartlett, W.A., Hewison, M., Stewart, P.M., 2002. Osteoblastic 11 β-hydroxysteroid dehydrogenase type 1 activity increases with age and glucocorticoid exposure. J. Bone Miner. Res. 17 (6), 979–986. doi: 10.1359/jbmr.2002.17.6.979 .

Eurostat, 2008. Your key to European statistics . http://ec.europa.eu/eurostat/web/ conferences/conf-2008 . Online; accessed 26 Jun 2017

Evans, N., Chapman, M., Chappell, M., Godfrey, K., 2002. Identifiability of uncon- trolled nonlinear rational systems. Automatica 38 (10), 1799–1805. doi: 10.1016/ S0 0 05-1098(02)0 0 094-8 .

Feldman, D. , Langberg, M. , 2011. A unified framework for approximating and clus- tering data. In: Proceedings of 43rd Annual ACM Symposium on Theory of Com- puting (STOC 2011), pp. 569–578 .

Gathercole, L.L., Lavery, G.G., Morgan, S.A., Cooper, M.S., Sinclair, A.J., Tomlin- son, J.W., Stewart, P.M., 2013. 11ß-hydroxysteroid dehydrogenase 1: translational and therapeutic aspects. Endocr. Rev. 34 (4), 525–555. doi: 10.1210/er.2012-1050 .

Godfrey, K.R., DiStefano, J., 1987. Chapter 1 - identifiability of model parameters. In: Walter, E. (Ed.), Identifiability of Parametric models. Pergamon, pp. 1–20. doi: 10.1016/B978- 0- 08- 034929- 9.50 0 05-4 .

Hassan-Smith, Z.K., Morgan, S.A., Sherlock, M., Hughes, B., Taylor, A.E., Lavery, G.G., Tomlinson, J.W., Stewart, P.M., 2015. Gender specific differences in skeletal mus- cle 11 β-HSD1 expression across healthy ageing. J. Clin. Endocrinol.Metab. 100 (7), 2673–2681. doi: 10.1210/jc.2015-1516 .

Hassan-Smith, Z.K., Sherlock, M., Reulen, R.C., Arlt, W., Ayuk, J., Toogood, A .A ., Cooper, M.S., Johnson, A.P., Stewart, P.M., 2012. Outcome of Cushing’s disease following transsphenoidal surgery in a single center over 20 years. J. Clin. En- docrinol.Metab. 97 (4), 1194–1201. doi: 10.1210/jc.2011-2957 . PMID: 22278423

Hundertmark, S. , Buhler, H. , Rudolf, M. , Weitzel, H.K. , Ragosch, V. , 1997. Inhibition of 11 β-hydroxysteroid dehydrogenase activity enhances the antiproliferative effect of glucocorticosteroids on MCF-7 and ZR-75-1 breast cancer cells. J. Endocrinol. 155 (1), 171–180 .

Jacquez, J.A. , 1996. Compartmental analysis in biology and medicine. BioMedware, Ann Arbor, MI .

Ketchen, D.J., Shook, C.L., 1996. The application of cluster analysis in strategic man- agement research: an analysis and critique. Strateg. Manag. J. 17 (6), 441–458. doi: 10.1002/(SICI)1097-0266(199606)17:6  441::AID-SMJ819  3.0.CO;2-G . Kohonen, T. , 1982. Self-organization of topologically correct feature maps. Biol. Cy-

bern. 43, 59–69 .

Lloyd, S., 1982. Least squares quantization in PCM. IEEE Trans. Inf. Theor. 28 (2), 129–137 . Originally as an unpublished Bell laboratories Technical Note (1957) doi: 10.1109/TIT.1982.1056489 .

Markey, K., Uldall, M., Botfield, H., Cato, L., Miah, M., Hassan-Smith, G., Jensen, R., Gonzalez, A., Sinclair, A., 2016. Idiopathic intracranial hypertension, hormones, and 11ß-hydroxysteroid dehydrogenases. J. Pain Res. 9, 223-232. doi: 10.2147/ JPR.S80824 .

Morgan, S., Hassan-Smith, Z., Doig, C., Sherlock, M., Stewart, P., Lavery, G., 2016a. Glucocorticoids and 11ß-HSD1 are major regulators of intramyocellular protein metabolism. J. Endocrinol. 229 (3), 277–286. doi: 10.1530/JOE- 16- 0011 .

Morgan, S., Hassan-Smith, Z., Doig, C., Sherlock, M., Stewart, P., Lavery, G., 2016b. Mechanisms in endocrinology: tissue-specific activation of cortisol in Cushing’s syndrome. Eur. J. Endocrinol. 175 (2), R83–9. doi: 10.1530/EJE-15-1237 .

Morgan, S.A., McCabe, E.L., Gathercole, L.L., 2014. 11ß-HSD1 is the major regulator of the tissue-specific effects of circulating glucocorticoid excess. Proc. Natl. Acad. Sci. U.S.A. 111 (24), E24 82–E24 91. doi: 10.1073/pnas.1323681111 .

Newell-Price, J., Bertagna, X., Grossman, A.B., Nieman, L.K., 2006. Cushing’s syn- drome. Lancet 367 (9522), 1605–1617. doi: 10.1016/S0140- 6736(06)68699- 6 . Owen, J. , Fiedler-Kelly, J. , 2014. Introduction to Population Pharmacokinetic/

Pharmacodynamic Analysis with Nonlinear Mixed Effects Models. John Wiley & Sons .

Sager, J. , Yu, J. , Raguenau-Majlessi, I. , Isoherranen, N. , 2015. Physiologically based pharmacokinetic (PBPK) modeling and simulation approaches: a systematic re- view of published models, applications and model verification. Drug Metab. Dis- pos. 43, 1823–1837 .

Sakuma, K., Yamaguchi, A., 2013. Sarcopenic obesity and endocrinal adaptation with age. Int. J. Endocrinol. 2013 (204164). doi: 10.1155/2013/204164 .

Self, A. , Thomas, J. , Randall, C. , 2012. Measuring National Well-being, First Annual Report on Measuring National Well-being, 2012. Technical Report. Office for Na- tional Statistics .

van Staa, T., Leufkens, H., Abenhaim, L., Begaud, B., Zhang, B., Cooper, C., 20 0 0. Use of oral corticosteroids in the United Kingdom. QJM 93 (2), 105–111. doi: 10.1093/ qjmed/93.2.105 .

Tiganescu, A., Tahrani, A.A., Morgan, S.A., Otranto, M., Desmoulière, A., Abrahams, L., Hassan-Smith, Z., Walker, E.A., Rabbitt, E.H., Cooper, M.S., Amrein, K., Lav- ery, G.G., Stewart, P.M., 2013. 11ß-hydroxysteroid dehydrogenase blockade pre- vents age-induced skin structure and function defects. J. Clin. Invest. 123 (7), 3051–3060. doi: 10.1172/JCI64162 .

Tomlinson, J.W., Sherlock, M., Hughes, B., Hughes, S.V., Kilvington, F., Bartlett, W., Courtney, R., Rejto, P., Carley, W., Stewart, P.M., 2007. Inhibition of 11 β- hydroxysteroid dehydrogenase type 1 activity in vivo limits glucocorticoid expo- sure to human adipose tissue and decreases lipolysis. J. Clin. Endocrinol.Metab. 92 (3), 857–864. doi: 10.1210/jc.2006-2325 .