The Non-Linear Child: Ontogeny, Isoniazid Concentration, and NAT2 Genotype Modulate Enzyme Reaction Kinetics and Metabolism.

N-acetyltransferase 2 (NAT2) catalyzes the acetylation of isoniazid to N-acetylisoniazid. NAT2 polymorphism explains 88% of isoniazid clearance variability in adults. We examined the effects of clinical and genetic factors on Michaelis-Menten reaction kinetic constants of maximum velocity (Vmax) and affinity (Km) in children 0-10years old. We measured the rates of isoniazid elimination and N-acetylisoniazid production in the blood of 30 children. Since maturation effects could be non-linear, we utilized a pharmacometric approach and the artificial intelligence method, multivariate adaptive regression splines (MARS), to identify factors predicting NAT2 Vmax and Km by examining clinical, genetic, and laboratory factors in toto. Isoniazid concentration predicted both Vmax and Km and superseded the contribution of NAT2 genotype. Age non-linearly modified the NAT2 genotype contribution until maturation at ≥5.3years. Thus, enzyme efficiency was constrained by substrate concentration, genes, and age. Since MARS output is in the form of basis functions and equations, it allows multiscale systems modeling from the level of cellular chemical reactions to whole body physiological parameters, by automatic selection of significant predictors by the algorithm.

Introduction

Isoniazid is a first line antituberculosis agent recommended for treatment of children with tuberculosis (TB). Its conversion is catalyzed by arylamine N-acetyltransferase (NAT) isoenzyme 2 (NAT2) in phase II xenobiotic metabolism. NAT2 is central to detoxifying many environmental-, industrial-, and food-based arylamines and hydrazines. NAT2 catalyzes the reaction of acetyl-CoA and isoniazid (substrates) to produce CoA and N-acetylisoniazid. NAT2 single nucleotide polymorphisms (SNPs) identified in adults led to classification of individuals as slow and rapid acetylators. Specific SNPs could be associated with changes in enzyme catalytic activity, as defined by enzyme reaction kinetic constants in liver tissue (Blum et al., 1991, Fretland et al., 2001). The elimination of isoniazid from the body is a combination of NAT2 metabolism and direct elimination of parent compound, with parent isoniazid elimination in urine of up to 37%. However, in adult patients 88% of systemic clearance variability has been linked to NAT2 SNPs (Kinzig-Schippers et al., 2005). While the ontogeny of phase I enzymes is well established, the age-dependent maturation changes in phase II enzymes such as NAT2 are still poorly understood. This is important since isoniazid peak concentrations and 0–24 h area under the concentration-time curves (AUC0–24) are major determinants of cure rates, speed of sterilizing effect, and acquired drug-resistance, in TB patients (Chigutsa et al., 2015, Pasipanodya et al., 2013, Pasipanodya et al., 2012, Pasipanodya and Gumbo, 2011, Gumbo et al., 2014, Dheda et al., 2014). There could be an age-dependent difference in systemic clearance of isoniazid in children, however, the effect of maturation, or other pediatric factors, on NAT2 enzyme kinetics has hitherto not been investigated (Jeena et al., 2011, Rey et al., 2001). Here, we identified the rate of production of N-acetylisoniazid and the rate of elimination of isoniazid in children enrolled in the prospective PHATISA study (Hiruy et al., 2015) and investigated the possible effects of NAT2 SNPs, measures of body-size, age, and nutrition status, in altering NAT2 reaction kinetics and isoniazid elimination rate constants.

Human growth and physiological maturation are non-linear processes, with discordant changes that have direct effect on drug pharmacokinetics (Kearns et al., 2003). Indeed, in recent years it has become evident that even in adults the relationships between different pharmacokinetic parameters and covariates such as age and weight are not only non-linear, but have high order interactions, with data zones of discontinuity in the relationships (Hall et al., 2011, Hall et al., 2013, Jain et al., 2013). Thus, standard statistical analyses could be limited on several fronts, including reliance on linear analyses and the restrictive normality, collinearity and sampling assumptions needed for valid inference. In 2012, Kiser et al. used standard statistical and pharmacokinetic approaches and identified that NAT2 genotype and child's age were associated with isoniazid clearance in South African infants. With the slow acetylator genotype, weight adjusted clearance was unchanged from week 0 to week 12 but increased 28% in week 84; while for rapid acetylators it increased 13% from week 0 to week 12 but then decreased by 23% in week 84 (Kiser et al., 2012). In other words, the changes were non-linear. Isoniazid elimination is from a combination of NAT2 metabolism and loss of parent compound in urine, which could account for the complex picture observed. Here, we used non-linear science to identify the role of maturation, weight, and NAT2 SNPs on xenobiotic metabolism in children (Campbell, 1987, Campbell et al., 1985, Dokoumetzidis et al., 2002). We applied artificial intelligence (AI)-based non-linear analyses to investigate NAT2 maximum velocity (Vmax) and affinity (Km) during the first decade of life in children with TB. We were interested in determining relationships that could be used to scale from the level of enzyme-based chemical reactions to the whole organism (patient) in the clinic, to allow for direct translation.

Machine learning is a branch of AI-based methods that we have applied to examine pathophysiological and clinical parameters as well as therapeutic outcomes in the past, whereby it outperformed standard statistical inference with datasets of as little as 28 patients (Chigutsa et al., 2015, Gumbo et al., 2014, Modongo et al., 2015, Modongo et al., 2016, Jain et al., 2013, Pasipanodya et al., 2013, Pasipanodya et al., 2015). One method, multivariate adaptive regression splines (MARS), was introduced by Friedman (1991). MARS has been utilized by several groups for a variety of purposes, with both large and small clinical sample sizes for a larger number of predictors (Brasier et al., 2012, George and Chang, 2014, Lee et al., 2014, Lin et al., 2008, Liu et al., 2015, Ju et al., 2014). Studies included examination of concentrations of 8 proteins as potential predictors in 30 patients for helicobacter infection, RNA-sequencing data to differentiate genes with low and high expression in the entire human genome from sample sizes as low as four per group, and an examination of 419 plasma protein concentrations as potential predictors of dengue fever in 55 patients. MARS uses an optimization procedure that employs basis functions (BF) to examine candidates for both main and interaction effects in a piecewise manner. BFs are flexible mathematical linear spline functions that describe relationship between variables in short segments of dataspace, which allows fit of both nonlinear and linear relationships as well as their interactions simultaneously. A backward elimination step is used to select a set of BFs that minimizes mean-square errors, after which multifold cross validation procedure is employed for model validation. The general cross-validation (GCV) step gives a measure of how well the model will perform on an independent dataset. In AI, an interaction is encountered when one predictor modifies another predictor's effect on the outcome (say Vmax). In the same model, MARS identifies inflection thresholds (termed “hinges”), where relationship between predictor and target changes or even ceases to exist within specific and well-defined data-space. This essentially allows the slope of the relationship between target and predictor to be different for different data ranges of the same predictor variable, and allows for discontinuity of relationships in some data spaces. Finally, the variable importance score is used for ranking and selecting predictors, especially important when the number of predictors is relatively large compared to the number of observations (as was the case in our study) and when many of the predictors could be co-linear (Breiman, 2001). Here, we applied MARS to data obtained from an experimental clinical design to identify predictors of NAT2 Vmax and Km in South African children with TB.

The philosophy and purpose of AI methods such as MARS differ from standard statistical inferences (Breiman, 2001, Hastie et al., 2009). The underlying scientific philosophy in standard statistical inference is rejection of a pre-specified null hypothesis. Comparison of distributions and the measures of central tendency (or hazard rates) between two or more groups are often the main way to test hypotheses, with tests of significance designed to make sure that differences observed are not due to chance. The purpose of MARS is not for hypothesis testing, but rather for hypothesis generation. Since these algorithms are distribution free, measures of central tendency are not the issue. Moreover, MARS algorithm uses non-parametric regression: in other words, predictors do not take a pre-specified form but are derived from the data. The main objective of MARS is prediction and pattern recognition, based on present data. Thus, there is no test of significance. Rather, the algorithms determine how likely the prediction is likely to be true with different datasets. During the GCV-step, the algorithms randomly split the dataset into different sized smaller sets and calculate a GCV score, which tells us how well predictors will perform on an entirely new data-set. When predictors and their interactions are identified with these approaches, they become hypotheses to be tested using standard statistical inference.

Methods

PHATISA Study Design

The study design, regulatory approvals, and the summary pharmacokinetic parameter results of the PHATISA study have been published in the past (Hiruy et al., 2015). Clinical and demographic factors were recorded. Malnutrition was measured using Z scores for height-for-age, weight-for-age, and weight-for-height, and mid-upper arm circumference (MUAC), since it has been proposed that isoniazid elimination could be affected by malnutrition. The nutrition measures such as MUAC and Z-scores were interpreted based on WHO criteria (http://whqlibdoc.who.int/hq/1997/WHO_NUT_97.4.pdf). The children had blood draws after directly observed therapy between the fourth and twelfth day after initiation of anti-TB therapy. In order maximize experimental design and limit the number of samples drawn, and to reduce the number of children under study, times for blood draw were specifically identified using optimal sampling theory (Reed, 1999, Wang and Endrenyi, 1992). We drew blood for concentration of both isoniazid and its metabolite, N-acetylisoniazid at 0, 0.42, 1.76, 3.37, 10.31 and 24 h after an isoniazid dose. Blood was processed for extraction of isoniazid and N-acetylisoniazid, as well as for DNA extraction.

Isoniazid and N-acetylisoniazid Concentration Measurement Assays

N-acetylisoniazid was synthesized from isoniazid standard using acetic anhydride and a previously published method. Following sample preparation using an acetonitrile-based protein precipitation process, measurement of drug and metabolite concentrations was carried out by liquid chromatography-tandem mass spectrometry on an AB Sciex 5500 Q-trap mass spectrometer coupled to an Agilent 1200 UPLC system. An Exterra 2.1 mm × 50 mm C18 column and a mobile phase of water with 0.1% formic acid at a flow rate of 200 μL/min was used for separation. The mass spectrometer was operated in the positive ion mode and the following transitions monitored: Isoniazid: 138.1 → 51.9 and 138.1 → 66, N-acetyl isoniazid: 180 → 66 and 180 → 78.6 and 6-aminonicotinic acid (internal control): 138.7 → 51 Analyst® 1.5 software version 1.5.1 was used for quantitation.

NAT2 Genotyping

DNA was extracted from buffy coats using the Qiagen Flexigene DNA Kit. The entire NAT2 gene was amplified by PCR using the sense primer (5′ GAC ATT GAA GCA TAT TTT GAA AG 3′) and the antisense primer (5′ GAT GAA AGT ATT TGA TGT TTA GG 3′) primer pair (Hickman and Sim, 1991). Amplification conditions were as follows: an initial denaturation of 5 min at 95 °C, followed by 35 cycles of denaturation at 95 °C for 30 s; annealing at 55 °C for 30 s; extension at 72 °C for 60 s, with a final extension of 5 min at 72 °C. Sanger sequencing (performed by Inqaba Biotechnical Industries) was used to determine NAT2 genotypes. NAT2 alleles were characterized by polymorphisms located at seven positions 191, 282, 341, 481, 590, 803, and 857, as described by the Database of arylamine N-acetyltransferases (NATs) (http://nat.mbg.duth.gr; McDonagh et al., 2014, Sabbagh et al., 2011). SNPs at these 7 positions have a genotype/acetylation status prediction of accuracy of 98.4% (Cascorbi et al., 1995, Hein and Doll, 2012).

Identification of Enzyme Kinetics and Pharmacokinetic Parameter Estimates

We implemented two types of models from the ADAPT software of D'Argenio et al. (2009). The parameter estimates for each child were estimated using the maximum-likelihood solution via the expectation-maximization algorithm (MLEM). The changes in concentration of isoniazid (INH) and N-acetylisoniazid (ANIH) with time were utilized in two types of pharmacometric models. The first was a pharmacokinetic model for the plasma concentrations of parent drug and its metabolite, with no assumptions of NAT2 saturability. The isoniazid elimination rate (KINH) as well as the fraction of isoniazid metabolized (fINH) and the volume of distribution of metabolite VAINH were used to calculate the elimination ratio VAINH/fINH. In order to identify NAT2 Vmax (mg/h) and Km (mg/L) for isoniazid metabolism in each child, we co-modeled the concentration versus time of isoniazid and N-acetylisoniazid using a 2-compartment model, chosen based on our prior work with this drug (Hiruy et al., 2015, Pasipanodya et al., 2013). A standard two-stage estimation method was used to generate initial Km and Vmax values and other pharmacokinetic parameter estimates. The parameters identified were then used in subroutine POPINIT of ADAPT, and enzyme reaction kinetic constants and pharmacokinetic parameter estimates for each child were then identified using MLEM.

Frequentist Statistical Analysis

We employed the D'Agostino and Pearson omnibus normality test and a p-value > 0.05 for normality tests. Comparison between proportions was made using the Fischer's exact or χ-square tests. Comparison between medians or means made via the Kruskal-Wallis test. All statistical testing other than machine learning methods were performed using STATA version 13 software (STATA Software, College Station, Texas).

Machine Learning

MARS uses BFs to examine candidates for both main and interaction effects, which gives MARS flexible and adaptive capabilities to fit nonlinear and linear relationships and interaction components simultaneously. Over fitting procedures were used to grow large models with up to 15 BFs, which were then pruned back using GCV function during the backward pass. Each variable was assigned a measure of predictive importance by MARS, entailing both marginal and interaction effects involving this variable. We examined for predictors of four target phenotypic responses: Vmax, Km, KINH and VAINH/fINH. The potential predictors entered for each model for each of the four responses were each child's NAT2 genotype, age in days, INH dose, observed peak concentration, HIV status, gender, height, weight, body mass index (BMI) in kg/m2, Z-scores for both height and weight, as well as MUAC. The coefficient of determination (R2) values, variable importance measures and the MARS equation for selected models are reported. All machine learning modeling was performed on the SPM Salford Predictive Modeler® software suite version 7.0 (SALFORD Software, San Diego, California). Except for setting maximum allowable basis functions to 15, all other settings were set to default for both the exploratory legacy mode and in cross-validation.

Results

We enrolled children ≤ 10 years old who presented with TB to King Edward VIII Hospital in Durban, South Africa, between May 2012 and March 2013. We enrolled 19 boys and 11 girls, whose clinical characteristics are shown in Fig. 1. The figure shows that most children were one year-old or less. All children had chest X rays performed, 19 (63%) had consolidation and 4 (13%) had cavitation, thus most children had pulmonary TB. Fig. 2 shows Z-scores and demonstrates that up to 30% of children were malnourished, though only one was severely malnourished (Z-score standard deviation < − 3).

Fig. 1

Demographic and dosing characteristics in 30 children. The p-values shown in each graph are for the D'Agostino and Pearson omnibus normality test. A p-value > 0.05 is significant, while < 0.05 is non-significant for normal distribution. (A). Age was the only parameter not normally distributed and showed a marked skew toward the younger age groups. This in turn dictated the weight (B) and height (C) at lower values of these, given the link between age and these parameters. (D). The isoniazid mean ± standard deviation dose was 11.5 ± 4.9 mg/kg; consistent with the new World Health Organization guidelines (2010).

Fig. 2

Z-scores and mid-upper arm circumference as indicators of nutrition in 30 children (A). Z-score for height is shown, with 4/30 (13%) children designated as malnourished by that measure alone. (B). Z-score for weight means that 7/30 (23%) children were malnourished by this measure. (C). Height for weight score suggested that 9 (30%) children could classified as malnourished. (D). The mid-upper circumference distribution suggests that 3/30 (10%) of children were malnourished.

Isoniazid acetylation follows a ping-pong Bi-Bi mechanism described by the standard Michaelis-Menten relationship for enzyme rate reaction rate (v) and substrate concentration [S] (Michaelis and Menten, 1913, Weber and Cohen, 1968):

We generated Km and Vmax values for each child in ADAPT based on this relationship, with estimates shown in Fig. 3A–B. The mean population estimate for Km was 1.69 mg/L, which translates to 0.012 μM, while that of Vmax was 30.11 mg/h, which translates to 0.06 μM/s. Fig. 3A-B show that Km was tightly distributed around the mean, but there was greater variability in the Vmax. The isoniazid elimination rate (KINH) is shown in Fig. 3C. Fig. 3D shows the elimination ratio VAINH/fINH, calculated using the fraction of isoniazid metabolized (fINH), and the volume of distribution of metabolite VAINH.

Fig. 3

N-acetyltransferase enzyme kinetic and isoniazid elimination constants in children with tuberculosis. Mean values with (±) standard deviations for all 30 children. (A). Maximum velocity (mg/h) value distribution (B). Affinity (mg/L). (C). Isoniazid elimination rate constant (per hour). (D). Elimination ratio distribution.

NAT2 sequencing results, including the base pairs and SNP calls for each allele, are shown in Extended Dataset 1. We identified 17 different alleles, including four reference high-activity NAT2*4 alleles. The most frequent genotypes were NAT2*5B/NAT2*5B encountered in 6/30 (20%) children, and NAT2*5C/NAT2*12C encountered in 6/30 (20%) children; the former genotype is commonly reported in literature to confer slow acetylator status, while the latter is predicted to confer the intermediate acetylator status (McDonagh et al., 2014). On the other hand, homozygosity for the rapid acetylator NAT2*4 reference allele (NAT2*4/NAT2*4 genotype) was encountered in 2/30 (7%) of enrolled children (Table 1).

Table 1

List of NAT2 genotypes and different phenotypes in 30 tuberculosis children.

NAT2 Genotype	Inferred acetylation status	Number of children	Isoniazid pharmacokinetic parameter estimates	NAT2 enzyme function
			Isoniazid elimination rate (KINH)	Maximum velocity (Vmax)	Affinity (Km)
*4/*4	Rapid	2	10.306 (0.006)	30.460 (0.335)	1.678 (0.010)
*4/*5C	Intermediate	1	10.307	28.797	1.716
*4/*14E	Intermediate	2	10.306 (0.003)	29.919 (1.733)	1.699 (0.049)
*5B/*5B	Slow	6	10.311 (0.016)	29.563 (3.679)	1.717 (0.067)
*5B/*5C	Slow	2	10.308 (0.018)	30.349 (0.608)	1.687 (0.018)
*5C/*5K	Slow	1	10.306	29.640	1.699
*5C/*6N	Slow	2	10.308 (0.006)	30.802 (1.578)	1.677 (0.043)
*5C/*11A	Intermediate	1	10.320	30.198	1.687
*5C/*12A	Intermediate	1	10.308	29.823	1.700
*5C/*12C	Intermediate	6	10.311 (0.012)	30.233 (1.395)	1.683
*5C/*14E	Slow	1	10.309	30.163	1.690
*5D/*11A	Intermediate	1	10.311	30.864	1.670
*6A/*12A	Intermediate	1	10.301	30.864	1.670
*6B/*6B	Slow	1	10.304	28.231	1.758
*12A/*12B	Rapid	1	10.292	31.690	1.652
*12A/*13A	Rapid	1	10.317	31.881	1.646

As a first step, we employed standard frequentist statistics to compare the distributions of Km, Vmax, and elimination rate between different inferred genotypes and predicted acetylation status. Table 2 shows no significant differences in mean values. Compared the distributions of the same enzyme function parameters, as well as measures of body size, age, and dose by inferred acetylator status, shown in Table 2. The only statistically significant finding was that the isoniazid dose administered was lower among rapid acetylators.

Table 2

Distribution of demographic characteristics and measures of NAT2 enzyme function between inferred acetylation status groups.

	Inferred acetylation status
	Rapid (n = 3)	Intermediate (n = 13)	Slow (n = 14)	P-valuea
Enzyme function
Isoniazid elimination	10.301 (0.010)	10.310 (0.010)	10.308 (0.012)	0.2945
Maximum velocity (Vmax)	30.870 (0.749)	30.206 (1.234)	29.853 (2.422)	0.273
Enzyme affinity (Km)	1.669 (0.017)	1.686 (0.028)	1.704 (0.049)	0.2122
Demographic characteristics
Age-in-days	1131.5 (867.3)	1151.4 (942.5)	1064.76 (1094.61)	0.7297
Weight (kg)	12.67 (5.13)	12.83 (5.25)	12.06 (6.01)	0.8655
Height (cm)	90 (18.35)	90.54 (16.68)	85.14 (20.85)	0.6456
Body mass index (kg/m2)	15.20 (2.29)	15.25 (2.65)	15.95 (4.04)	0.9462
Isoniazid dose (mg/kg)	6.29 (1.98)	15.80 (4.40)	11.39 (5.22)	0.099
Isoniazid total daily dose (mg)	73.33 (15.28)	147.69 (41.06)	130.71 (102.37)	0.017

Based on K-Wallis test.

Next, we used MARS algorithm in order to identify predictors of Vmax. Table 3 shows that dose (BF1), genotype (BF5), isoniazid peak concentration (BF7), and age (BF10) were identified as the best predictors of Vmax. The table illustrates the nature of MARS output, which is the form of basis functions (BF). BF1 is a simple hinge function “max (0, Dose–190)” which means that the value of the expression is zero for all doses (mg) up to 190 (i.e., for all doses satisfying “Dose minus 190” less than zero), as illustrated in Fig. 4A. The dose of 190 mg is at the hinge. The variable importance score of 100% means that “Dose” is the most important variable for Vmax. Isoniazid peak concentration was also an important predictor (BF7), with a variable importance score of 90%, very close to that of “Dose”; thus both dose and peak concentration likely reflect the same thing, which is effect of isoniazid concentration on Vmax. BF7 means that isoniazid concentration had no effect on Vmax until the peak concentration achieved the threshold value of 8.83 mg/L, at the hinge. Table 3 shows that the peak concentration (BF7) interacted with BF5 (genotype). However, the BF5 interaction with BF7 only applied to NAT2 genotypes selected as subset 1, which means that there is staggered or stratified effect so that relationship applied to some genotypes and not others. Specifically, BF5 comprised NAT2 genotypes that were grouped by MARS, which had a variable importance score of 90%, similar to that of isoniazid concentration. Of note, for MARS analyses we did not pre-categorize different genotypes by acetylator status (slow or rapid acetylator), since predictors do not take a pre-specified form but are derived from the data. Table 3 further shows that, in addition, the child's age also had an effect on Vmax as well, as reflected by BF10; the impact however was less than half that of either genotype or concentration. Fig. 4B shows the singular contribution of age; the effect of age starts at zero but reaches a maximum at 1935 days (about 5.3 years) when it reaches a plateau, and remains flat into the next decade. Overall, the relation between Vmax and these four predictors was given by the equation:

Table 3

Factors predictive of maximal velocity (Vmax) on MARS analysis.

	Function	Units	Variable importance score (%)
BF1	max (0, Dose − 190)	Dose in mg	100
BF5	Genotype subset 1	Allele 1/Allele 2	90
BF7	max (0, Peak isoniazid concentration − 8.83)*BF5	mg/L	90
BF10	max (0, 1934.5-Age)	Age in days	42

NAT2 genotype subset #1: 5D/*11A, 5C/12C, 5C/*14E, *4/*4, *4/*5, *4/*5C, *5B/*5B, *5C/*12A, *5C/*5K, *6B/*6B.

Fig. 4

Contribution of dose and age on maximum velocity. (A). Effect of dose on Vmax, shows lack of effect of dose below 190 mg, after which dose has a positive relationship with Vmax. In this case, dose is a surrogate of drug concentration. (B). The figure summarizes the hinge function “max (0, 1934.5-Age)”, a perfect example of the flexibility of MARS in modeling change in relationship between predictor and outcome at a hinge, and thus non-linear. The figure shows that there is a large effect of age until 5.3 years when the contribution gets to a plateau, suggesting the point of full maturation.

Thus, age (maturation) affects speed of isoniazid acetylation by modifying effect of genotype, a non-linear interaction, however after 5.3 years the maturation is complete.

Table 4 shows the BFs for NAT2 affinity for isoniazid, or Km. A child's height was the primary predictor and had the highest variable importance score. Height however could be substituted by age, given the tight link between height and age (Pearson r = 0.917; p < 0.001), a tightly linked relationship frequently used in the clinic. Thus age was the main predictor of Km as well; in fact, the coefficient of determination (R2) value of age alone was 31%. However, MARS also identified that isoniazid peak concentration and genotype as important predictors (Table 4). Combined, these factors explained 80% of the variance in NAT2 enzyme affinity, as given by equation for Km and its predictors:

Table 4

Factors predictive of enzyme affinity (Km) on MARS analysis.

	Function	Units	Variable importance score (%)
BF1	max (0, height-54)	centimeters	100
BF2	Max (0, isoniazid peak concentration-11.7)	mg/L	53
BF4	Genotypes subset 1*BF1	Allele 1/Allele 2	26
BF7	Genotypes in subset 2	Allele 1/Allele 2	26

NAT2 genotype subset 1: *5D/*11A, *5C/*12C, *4/*14, *4/*5, *4/*5C, *5B/*5B, *5B/*5C, *5C/*5K, *5C/*6N.

NAT2 genotype subset 2: *12A/*13A, *6A/*12A, *4/*4.

Table 4 shows that BF4 is a special type of a basis function, described by an interaction of BF1 (age or height) and genotype. In other words, age constrained the effect of genotype on enzyme affinity.

Table 5 shows the BFs for isoniazid elimination constant (KINH). The elimination rate constant is a composite measure reflecting both acetylation and elimination of parent compound via loss in urine. The main predictors were isoniazid dose and NAT2 genotype. In this case, however BF2 was an interaction based on the modification of effect of genotype on the effect of isoniazid dose (i.e., isoniazid concentration), an example of higher order interaction among predictors. This relationship is shown in Eq. (4):

Table 5

Factors predictive of isoniazid elimination rate constant.

	Function	Units	Variable importance score (%)
BF1	max (0, Dose − 40)	mg	100
BF2	Genotypes subset 1*BF1	Allele 1/Allele 2	19

NAT2 genotype subset 1: *5C/*12A, *12A/*13A, *5C/*12/, *5C/*14E, *5C/*12A, *5C/*6N, *6B/*6B.

NAT2 genotype subset 2: *5C/*11A, *12A/*12B, *6A/*12A, *4/*4, *4/*5, *4/*5C, *5B/*5B, *5B/*5C, *5C/*5K.

This was confirmed by standard regression methods, which showed that for each unit increase in isoniazid dose above 40 mg, KINH increased by an additional 2.82 × 10− 4 per hour for subset 1 versus 1.05 × 10− 4 per hour for subset 2 genotypes, a 2.69-fold change in slope.

As regards to the ratio VAINH/fAINH, the R2 for the relationship was only 0.55 in the best model. We thus considered that our potential predictors such as clinical and demographic factors were not as explanatory of this parameter, and thus no further explorations were made.

Discussion

In pharmacometrics, primary pharmacokinetic parameters are the observed values such as elimination rate constant (KINH), volume of distribution (both central and peripheral), and in the current study also Km and Vmax. These observations quantify rates and extents of underlying physicochemical processes. In other words, these parameters, including the enzyme reaction constants, are the biological phenotypes. On the other hand, in pharmacogenomics, inferred acetylation status is commonly termed the “acetylation phenotype.” Acetylation phenotype is a grouping of drug clearances (into 3 overlapping Gaussian curves). Moreover, commonly in clinical practice, acetylation status is assigned using 2-h post-dose isoniazid concentrations. However, the isoniazid elimination rate itself is a product of at least two physicochemical processes, acetylation and the loss of unmetabolised (parent) compound in urine. “Phenotyping” methods such as the 2-h concentration or concentration-time profiles or AUC0–24 are even further removed from the primary process of N-acetylation: they are products of several primary biological processes, namely absorption rate constant (and closely related time to maximum concentration), volume of distribution, and elimination rate constants. Only the last of these is related to acetylation, and only partially. The primary biological process of acetylation itself (i.e., the phenotype), sensu strictiore, is defined by Km and Vmax values, which are the “observed” values based on calculations using the Michaelis-Menten Eq. (1). These values are thus not inferred from genotypes, but are the “observed”. Thus, effect of genotype should be inferred based on these phenotypes, and not vice-versa. Here, we simply asked: what are the physiological and clinical factors that could affect these quantitative measures of phenotype, if at all, and which factors affect them the most? In terms of the genes, what groupings of NAT2 genes categorize the phenotypes of isoniazid NAT2 affinity and maximum velocity in children, during periods of highest growth? MARS, in which predictors do not take a pre-specified form or groupings (e.g. inferred acetylation status), is excellent at such an agnostic function, and utilizes only the presented data for its own categorization.

To answer these questions, we identified the Km and Vmax values for each child from clinical data, based on a pharmacometric approach. The mean population estimates of 0.012 μM for Km for isoniazid. In biochemical assays for measurement of Km that examined isoniazid acetylation by human NAT2 from 2-h post-mortem liver homogenates of a rapid acetylator in the laboratory, the Km values were 0.018 ± 0.004 μM (Weber and Cohen, 1968). This means that the pharmacometric approach of using the rate of production of N-acetylisoniazid and disappearance of isoniazid is a robust enough approach to calculate these values for each child, and gives values similar to those from enzyme kinetic experiments of liver biopsy specimens. In other words, the measurement of the phenotype was likely accurate. In the case of the NAT2 alleles we identified in the children, genotypes such as NAT2∗ 4, NAT2∗ 12A, NAT2∗ 12B, NAT2∗ 12C, and NAT2∗ 13 alleles, and the predominance of the NAT2*5 allele have been shown to be common in Black South African adults (Dandara et al., 2003, Loktionov et al., 2002). Thus, the distribution of the common alleles in our small pediatric cohort mirrors that in adults. Application of AI methods to these parameter estimates and genotypes identified several important patterns between the predictors and enzyme reaction kinetic constants. Our approach is tractable and follows steps that are easy to apply to any pediatric drug whose metabolite is known for each child, and for a population of children. Thus, our approach could be used to examine the effects of maturation and drug doses on both compartmental pharmacokinetic parameters and Michaelis-Menten relationships for other phase II metabolism reactions, an unexplored pediatric space (Kearns et al., 2003).

Second, a current problem in systems pharmacology is the difficulty in scaling mathematical models from the level of chemical reactions inside the cell to the level of the individual patient, and to populations of patients. The MARS equations output could allow this to be accomplished, starting with reaction kinetics at the level of a Michaelis-Menten reaction to whole body physiological parameters in the child (isoniazid elimination rate), based on relevant clinical and demographic predictors automatically selected by the algorithm. Moreover, MARS handles both linear and non-linear interactions simultaneously, allowing for construction of non-linear models. Furthermore, MARS is assumption and distribution free and is flexible enough to identify interactions without a prior falsifiable hypothesis to link reactions constants to a predictor, as is the case of isoniazid concentration and Km. The findings led to new hypotheses being generated, which can now be tested in separate studies using standard statistical approaches. Thus, the mathematical relationships are identified first, and then physiological meanings can be investigated further to improve the precision of the estimates for better reproducibility.

Third, we found that the relationships between predictors such as age, dose or isoniazid concentration, and NAT2 enzyme kinetics were non-linear. In other words, several clinical and physiological factors affected speed of isoniazid acetylation and affinity, and interacted with each other. Of major importance was the effect of age on both speed of acetylation and enzyme affinity. Age played a critical role in modulating these enzyme activities, and age's contribution changed with child's age, with the maximal adult type of activity encountered after 5.3 years old. This time period likely reflects the time to maturation of the NAT2 enzyme. On the other hand, the effect of dose and drug concentration always superseded that of NAT2 genotype. The effect of dose and drug concentration was persistent regardless of the enzyme reaction constant examined, and was often the primary predictor, except in the case of Km. The physiological and pharmacological basis for this remain to be worked out, but the fact that the hinge was always encountered beyond a certain concentration suggests a real concentration-dependent effect, consistent with many pharmacological processes. The GCV step identified R2 ≥ 80% on post-test sets, which means there is a high likelihood of the same observations with future data-sets, so that the findings are likely valid. However, clearly nothing in our current understanding of Km and Vmax can explain how the pharmacological and clinical factors we identified as predictors could affect both enzyme affinity and maximum velocity for isoniazid acetylation. However, the historical disaster of gray baby syndrome in children treated with the drug chloramphenicol suggests that maturation has a major effect on enzyme function (Glazer et al., 1980, Mulhall et al., 1983, Sutherland, 1959). The functional deficiency in the UDP-glucuronyl transferase enzyme system encountered in the first few weeks of life, especially in premature neonates, led to cardiovascular collapse when high concentrations of the antibiotic were achieved in the children (Glazer et al., 1980, Mulhall et al., 1983, Sutherland, 1959). However, we reiterate that, the philosophy behind machine learning methods such as MARS is prediction and pattern recognition, and not identification of causal pathways. Thus, the next step is to test the effect of age and isoniazid concentration on NAT2 Km and Vmax using standard hypothesis testing approaches, and if confirmed then a new understanding of NAT2 catalyzed reaction kinetics should be sought.

Fourth, our findings could impact the most widely used pharmacokinetic models of maturation as commonly recommended by regulatory authorities (European Medicines Agency, 2007). Several methods are used to examine the effect of maturation on pharmacokinetic parameters, with the two most common being allometric scaling based on the ¾ fractal geometry law and physiologically-based pharmacokinetic modeling (Hope et al., 2007). These methods use specific drug clearance rates derived in adults to estimate drug clearances in children, with ontogeny effects scaled based on dry tissue weight of the organ clearing the drug. The central assumption is that xenobiotic metabolism enzyme reaction rates are the same between adults and children, and that what varies is the quantity of the enzymes proportional to organ weight and size. It is recognized in some cases, as for example sirolimus, that there could be different fetal isoform enzymes that get replaced by others after birth, thereby leading to altered clearance (Emoto et al., 2015); nevertheless the reaction kinetic constants of each of the different isoforms are assumed fixed. Moreover, the development of clearance with maturation is described by a single mathematical function. Our data suggests that non-linear functions and analyses, with no a priori assumptions, could be more accurate.

Finally, we found no quantitative relationship between several measures of malnutrition such as Z-scores and isoniazid metabolism. Indeed, no single model included Z-scores and MUAC as a predictor. However, most children had moderate malnutrition, and it could be that effects are encountered only with severe malnutrition. Thus, this question remains unresolved.

Our study has several limitations. The first is the small number of children enrolled. However, we deliberately optimized experimental design to identify information rich pharmacokinetic sampling so that the smaller numbers were deliberate. Moreover, MARS has successfully identified predictors in smaller populations than these ones in several studies by us and by others. Second, our findings lack physiological explanations. These will be subject of future studies. Third, the intent of AI approaches we used is not hypothesis testing, but rather pattern recognition. Thus, future studies in which our findings are tested as a falsifiable hypothesis are recommended.

The following are the supplementary data related to this article.

Author Contributions

Conception and design: T. Gumbo, B.M., J. Pasipanodya, P. Jeena, W. Bishai.

Development of methodology: Z. Rogers, H. Hiruy, J. G. Pasipanodya, J. Adamson, P. Jeena, W. Bishai, T. Gumbo

Acquisition of data: Z. Rogers, H. Hiruy, C. Mbowane, J. Adamson, L. Ngotho, F. Karim, P. Jeena, W. Bishai

Analysis and interpretation of data: Z. Rogers, J. Pasipanodya, P. Jeena, W. Bishai, T. Gumbo

Writing, review and/or revision of the manuscript: Z. Rogers, H. Hiruy, J. Pasipanodya, P. Jeena, W. Bishai, T. Gumbo

Study supervision: P. Jeena, W. Bishai

Conflict of Interest

All authors declare no conflict of interest.