Keywords: CR Categories: I.3.3 [Picture/Image Generation]: Antialiasing; I.3.7 [Three-Dimensional Graphics and Realism]: Color, shading, shadowing, and texture anti-aliasing, filtering, procedural texturing, rendering, shading, solid noise, solid texturing, texturing

Solid noise is a fundamental tool in computer graphics. Surprisingly, no existing noise function supports both high-quality antialiasing and continuity across sharp edges. In this paper we show that a slicing approach is required to preserve continuity across sharp edges, and we present a new noise function that supports anisotropic filtering of sliced solid noise. This is made possible by individually filtering the slices of Gabor kernels, which requires the proper treatment of phase. This in turn leads to the introduction of the phase-augmented Gabor kernel and random-phase Gabor noise, our new noise function. We demonstrate that our new noise function supports both high-quality anti-aliasing and continuity across sharp edges, as well as anisotropy.

Filtering Solid Gabor Noise

Introduction

Solid texturing [START_REF] Perlin | An image synthesizer[END_REF][START_REF] Peachy | Solid texturing of complex surfaces[END_REF]] is a popular method for objects that are sculpted or carved out of a solid material (e.g., a marble statue). To avoid excessive storage requirements, solid or 3D textures are typically procedural, and are often based on procedural solid noise (e.g., Perlin Noise [START_REF] Perlin | An image synthesizer[END_REF]). To achieve high-quality rendering, solid textures must be properly anti-aliased, similarly to traditional textures.

In recent years, there has been renewed interest in the problem of anti-aliasing procedural textures. This has resulted in the introduction of several new noise functions that support filtering [START_REF] Hart | Antialiased parameterized solid texturing simplified for consumer-level hardware implementation[END_REF][START_REF] Cook | Wavelet noise[END_REF][START_REF] Goldberg | Anisotropic noise[END_REF];

1,2 e-mail: ares.lagae@cs.kuleuven.be, george.drettakis@inria.fr [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]]. However, despite these recent advances, no existing noise function supports both high-quality anti-aliasing and continuity across sharp edges. We illustrate this in Fig. 1. Please also refer to the videos in the supplemental material, which illustrate this more clearly. More specifically, Perlin noise [START_REF] Perlin | An image synthesizer[END_REF]] results in detail loss when filtered (Fig. 1(a.2)). Wavelet noise [START_REF] Cook | Wavelet noise[END_REF] integrates solid noise perpendicularly to the surface of the object, which introduces discontinuities at sharp edges (Fig. 1(b.1)), since the normal changes discontinuously. Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]] projects 3D points onto the surface of the object along the surface normal, which, similarly to wavelet noise, does not preserve continuity over sharp edges (Fig. 1(c.1)).

For an in-depth discussion and comparison of these noise functions, please refer to the recent survey of Lagae et al. [2010a].

In this paper we show that a slicing approach is required to preserve continuity across sharp edges, and we present a new noise function based on Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]] that supports anisotropic filtering of sliced solid noise. We individually filter the slices of Gabor kernels. This requires the proper treatment of the phase of the kernel. We therefore introduce a new Gabor kernel, the phaseaugmented Gabor kernel. This, in turn, leads to our new noise function, random-phase Gabor noise. We also discuss how our derivations result in a generalization of the Projection-Slice Theorem as used by [START_REF] Cook | Wavelet noise[END_REF]. We demonstrate that our new noise function supports both high-quality anti-aliasing and continuity across sharp edges, as well as anisotropy.

Motivation

Texturing the surface of an object with filtered solid noise is conceptually a two-step process. In the first step, the 3D solid noise is mapped onto the 2D surface. This is done either by slicing the solid noise, i.e., evaluating the noise at the surface of the object (Fig. 2(a)) [START_REF] Perlin | An image synthesizer[END_REF]], or by projecting the solid noise, i.e., integrating the noise perpendicularly to the surface of the object (Fig. 2(b)) [START_REF] Cook | Wavelet noise[END_REF]. In the second step, the noise on the 2D surface is filtered. This is typically done using frequency clamping [START_REF] Norton | Clamping: A method of antialiasing textured surfaces by bandwidth limiting in object space[END_REF]], a general approach, or using a method-specific approach.

Filtering the noise on the surface using frequency clamping works better if the power spectrum of the noise on the surface is bandpass. The relation between the power spectrum of the 3D solid noise and the 2D noise on the surface is given by the Fourier Slice Theorem [START_REF] Cook | Wavelet noise[END_REF], which states that projecting in the spatial domain corresponds to slicing in the frequency domain (Fig. 2(a,c)), and vice versa (Fig. 2(b,d)).

Perlin [1985] uses slicing in his noise, which is typically filtered using frequency clamping. However, filtering introduces an aliasing vs. detail loss trade-off (Fig. 1(a.2)), since the power spectrum of Perlin noise is not band-pass (Fig. 2(c)) [START_REF] Cook | Wavelet noise[END_REF].

Cook and DeRose [2005] use projection and frequency clamping in their wavelet noise. However, even though the power spectrum of the noise on the surface is band-pass (Fig. 2(d)), filtering does not fully solve the aliasing vs. detail loss trade-off. This is because frequency clamping works at the level of entire octaves at once (see supplemental material for the exact method we use), and because frequency clamping is limited to isotropic filtering (Fig. 1(b.2)). projecting solid noise introduces discontinuities in the noise on the surface at sharp edges (Fig. 1(b.1)), since the normal changes discontinuously (Fig. 2(b)). [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF] use projection and a filtering approach specific to Gabor noise in their surface Gabor noise. They project Gabor kernels onto the surface and filter the resulting 2D Gabor kernels analytically. This results in high-quality anisotropic filtering (Fig. 1(c.2)), since the filtering approach works at the level of individual Gabor kernels. However, similarly to wavelet noise, projection introduces discontinuities at sharp edges (Fig. 1(c.1)).

In this paper, we present solid random-phase Gabor noise, a new noise function based on Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]. Our new noise uses slicing, which preserves continuity at sharp edges (Fig. 1(d.1)). Since filtering is inherently a 2D operation (i.e., it depends on surface orientation) [START_REF] Heckbert | Fundamentals of Texture Mapping and Image Warping[END_REF]], we have to explicitly model the slicing of the 3D Gabor kernels to be able to filter the resulting 2D Gabor kernels. This requires the introduction of a new Gabor kernel, the phase-augmented Gabor kernel, and a new Gabor noise, random-phase Gabor noise. Our new noise uses a filtering approach based on that of [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF], which supports high-quality anisotropic filtering (Fig. 1(d.2)).

In Sec. 3 we discuss slicing solid Gabor noise; the principal result is Eqn. 6, which shows how to slice a 3D phase-augmented Gabor kernel. In Sec. 4 we discuss filtering sliced solid Gabor noise; the principal result is Eqn. 10, which shows how to filter the resulting 2D phase-augmented Gabor kernel.

Slicing Solid Gabor Noise

In this section, we introduce a solid Gabor noise that supports slicing. We first show that the Gabor kernel of [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF] is not closed under slicing (Sec. 3.1). We consequently introduce the phase-augmented Gabor kernel, a new kernel that is closed under slicing (Sec. 3.2), and random-phase Gabor noise, a new noise using the phase-augmented Gabor kernel (Sec. 3.3). We finally show that solid random-phase Gabor noise supports slicing (Sec. 3.4). We also present a generalization of the Projection-Slice Theorem [START_REF] Cook | Wavelet noise[END_REF]

Slicing The Gabor Kernel

The Gabor kernel g of [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF] (Fig. 3

(a), image) is de- fined as g (x; a, ω) = e -πa 2 |x| 2 cos (2πω • x) , (1)
the multiplication of a Gaussian (Fig. 3(b), image) and a harmonic (Fig. 3(c), image), where R + is the bandwidth and ω ∈ R n is the frequency. Note that we wrote the kernel in vector notation and extended it to n dimensions. Slicing is denoted as the operator S subscripted by the hyperplane along which is sliced. We slice the kernel with a hyperplane Π (Fig. 3, solid green line), which results in

SΠ [g (x; a, ω)] = SΠ e -πa 2 |x| 2 SΠ [cos (2πx • ω)] , (2)
the product of a sliced Gaussian and a sliced harmonic. This is because S is a multiplicative linear operator (i.e., for every pair of functions f and g and scalar t,

S [f g] = S [f ] S [g] and S [tf ] = tS [f ]).
The hyperplane Π is specified by its unit normal n and its signed distance to the origin d (Fig. 3, dashed green line). Quantities in Π, such as the (n-1)-dimensional vector x ′ , are expressed relative to an arbitrary coordinate system, whose origin is the perpendicular projection of the origin of the n-dimensional space. We continue by slicing the n-dimensional Gaussian (Fig. 3(b), image), which results in

SΠ e -πa 2 |x| 2 x ′ = e -πa 2 d 2 e -πa 2 |x ′ | 2 , (3)
a weighted (n -1)-dimensional Gaussian (Fig. 3(b), red curve), where the weight e -πa 2 d 2 depends on the distance from the hyperplane to the center of the Gaussian. We also slice the n-dimensional harmonic (Fig. 3(c), image), which results in

SΠ [cos (2πx • ω)] x ′ = cos 2πproj Π ω • x ′ -2πdn • ω , (4)
a phase-shifted (n-1)-dimensional harmonic (Fig. 3(c), red curve), where proj Π ω is the perpendicular projection of ω onto Π, and the phase shift 2πdn • ω depends on both the distance from Π to the origin and the relative orientation of n and ω. Note that proj Π ω is the component of ω parallel to Π, and n • ω is the length of the component of ω perpendicular to Π. We now observe that, because of the phase shift (Fig. 3(a,c), φ s ), the product of the sliced Gaussian (Eqn. 3) and the sliced harmonic (Eqn. 4) cannot be expressed as a Gabor kernel (Eqn. 1). In other words, the Gabor kernel of [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF] is not closed under slicing.

Π φ s (a) a Π (b) ω Π φ s (c) φ Π φ s (d)

The Phase-Augmented Gabor Kernel

In order to overcome the problem with the phase shift, we now introduce the phase-augmented Gabor kernel, a new kernel that can accommodate the phase shift introduced by slicing.

Following the derivation of the previous section, we define the phase-augmented Gabor kernel g g (x; a, ω, φ) = e -πa 2 |x| 2 cos (2πx

• ω + φ) (5)
(Fig. 3(d), image), where φ ∈ [0, 2π) is the phase of the harmonic. Note that the Gabor kernel of [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF] is a special case of our kernel with zero phase (φ = 0), and that, in contrast to kernel of Lagae et al., Gabor filters used in image processing typically have a phase parameter as well. We now observe that slicing the n-dimensional phase-augmented Gabor kernel g results in

SΠ [g (x; a, ω, φ)] x ′ = w s g x ′ ; a, ω s , φ s , (6)
a weighted (n -1)-dimensional phase-augmented Gabor kernel (Fig. 3(d), red curve), where w s = e -πa 2 d 2 , ω s = proj Π ω and φ s = φ -2πdn • ω. Note that all quantities introduced by or affected by slicing are indicated with a superscript s . In other words, in contrast to the Gabor kernel of [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF], the phaseaugmented Gabor kernel is closed under slicing.

Random-Phase Gabor Noise

With the phase-augmented Gabor kernel, we can now introduce random-phase Gabor noise, a new noise that supports slicing.

The most important aspect of random-phase Gabor noise is how the phase parameter φ is handled. Since φ introduced only to ac-commodate slicing, it does not make sense to expose it as a noise parameter to the user, nor to fix it to a specific value, since it will be affected by subsequent slicing. We instead assign a random value to φ, i.e., the phases {φi} of the kernels are distributed according to a uniform distribution in the interval [0, 2π). The full motivation for this choice is provided in the next subsection. The phaseaugmented Gabor kernel with a random phase is called the randomphase Gabor kernel. The adjective random-phase is used when the phases are random, otherwise the adjective phase-augmented is used. We do not use the random weights {wi} used by [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF], since the noise has a zero mean due to the random phases. We thus define n-dimensional random-phase Gabor noise n as n (x; a, ω)

= i g (x -x i ; a, ω, φi) , (7)
where the random positions {x i } are distributed according to a ndimensional Poisson process with mean λ, and {φi} are the random phases. The statistical properties of a noise, i.e., its variance and power spectrum, determine its aspect, and are needed for several texturing operations, such as color mapping and filtering. We provide a complete set of analytical expressions for random-phase Gabor noise in the supplemental material. Note that our new noise results in a more elegant overall formulation in comparison with the noise of [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF] since many expressions are simplified.

Slicing Random-Phase Gabor Noise

In addition to the statistical properties of the 3D solid noise, texturing an object with solid noise also requires the statistical properties of the 2D noise on the surface, i.e., the sliced noise.

We show in the supplemental material that slicing an n-dimensional random-phase Gabor noise n with a hyperplane Π results in

SΠ [n (x; a, ω)] x ′ = √ w s n x ′ ; a, ω s , (8)
an (n -1)-dimensional random-phase Gabor noise, i.e., randomphase Gabor noise is closed under slicing. This means that the statistical properties of a sliced solid random-phase Gabor noise are obtained using the analytical expressions for the statistical properties of a 2D random-phase Gabor noise.

We were able to perform this derivation only when using random

The Slice-Projection Theorem for Stochastic Processes

By combining insight gained from Eqn. 8 with the definition of noise as a stochastic process [Lagae et al. 2010a, Sec. 2.2], we were able to generalize the Fourier Slice Theorem as used by [START_REF] Cook | Wavelet noise[END_REF]. More specifically, we generalize from a single noise instance and its Fourier transform to the stochastic process corresponding to the noise and its power spectrum.

Slice-Projection Theorem for Stochastic Processes. Let X be an n-dimensional stationary stochastic process with power spectrum SXX . Then the (n-1)-dimensional stochastic process Y = SΠ [X], which is the slice of X along the hyperplane Π, has the power spectrum SY Y = PΠ [SXX ], which is the projection of SXX onto Π.

For example, in Fig. 5, the 2D noise X (a, image) has the power spectrum SXX (b, image); therefore, the 1D noise Y (a, red curve), the slice of X along the line Π (a, green line), has the power spectrum SY Y (b, red curve).

Proof. The proof follows from three elements. (i) The fact that Y (the auto-correlation of Y ) is SΠ 0 [RXX ], which is the slice of RXX (the auto-correlation of X) along the hyperplane Π0 through the origin and parallel to Π. is because X is stationary, and independent of whether or not Π goes through the origin. (ii) The Wiener-Khinchin Theorem, which states that the auto-correlation and the power spectrum of a stochastic process are a Fourier transform [Papoulis and Pillai 2002, Ch. 9]. (iii) The Projection-Slice Theorem [Bracewell 1999, Ch. 16].

This theorem allows us to validate Eqn. 8, and to explain the aspect of a 2D slice of noise from the power spectrum of the 3D solid noise (e.g., in Fig. 4). Note that, although the Fourier Slice Theorem as used by [START_REF] Cook | Wavelet noise[END_REF] provides the basic intuition, this theorem is more general and more complete. For example, the proof of this theorem explains why the slicing plane Π does not have to go through the origin, which is not explained by the Fourier Slice Theorem alone. We believe this theorem might have applications beyond noise (e.g., 3D solid texture extrapolation from 2D exemplar slices).

Filtering Sliced Solid Gabor Noise

Now that we have a noise that supports slicing, we introduce a filtering approach for our new noise.

Filtering a sliced solid random-phase Gabor noise n with an filtering Gaussian G f is equivalent to filtering the sliced phase-augmented Gabor kernels,

SΠ [n (x; a, ω)] x ′ * G f x ′ ; 0, Σ f = i w s i g x ′ -x ′ i ; a, ω s , φ s i * G f x ′ ; 0, Σ f , (9)
where * denotes convolution, and G (x; µ, Σ) is a Gaussian with mean µ and variance matrix Σ. This is because S and * are linear operators, and because the phase-augmented Gabor kernel is closed under slicing (Eqn. 6). Note that n is a 3D function, and S [n], G f and g 2D functions on the plane Π. We thus derive an analytical expression for filtering a 2D phase-augmented Gabor kernel g with an anisotropic filtering Gaussian G f ,

g (x; a, ω, φ) * G f (x; 0, Σ f ) ≈ w f g x; a f , ω f , φ , ( 10
)
where

w f = cGF , a f ≈ 2π |ΣGF |, ω f = µGF ; ΣGF = Σ -1 G + Σ -1 F -1 , µGF = ΣGF Σ -1 G µG, cGF = cF G (µ G ; 0, ΣG + ΣF ); ΣF = 1 4π 2 Σ -1 f , cF = 1 2π |Σf | ; and ΣG = a 2 2π I, µ G = ω.
Note that all quantities introduced by or by filtering are indicated with a superscript f . Our filtering expression is a generalization of the one of Lagae et al. [2009, Eqn. 22, Eqn. 23]. Note that the approximate nature of Eqn. 10 (and Eqn. 11 below) is due to the approximation of the anisotropic filtered Gabor kernel with an isotropic kernel [Lagae et al. 2009, Eqn. 22, Eqn. 23 and Fig. 6], which works very well in practice.

We now have all the pieces in place to provide the expression for our new solid random-phase Gabor noise,

SΠ [n (x; a, ω)] x ′ * G f x ′ ; 0, Σ f ≈ i w s i w f i g x ′ -x ′ i ; a f , ω s,f , φ s i , (11)
where ω s,f is obtained by applying Eqn 6 and then Eqn 10 to ω. This means that our new noise is evaluated by summing a contribution for each 3D kernel close to the point of evaluation. This contribution is obtained by slicing the 3D kernel (Eqn. 6), filtering the resulting 2D kernel (Eqn. 10), and evaluating the filtered 2D kernel.

Implementation

The implementation of our new solid random-phase Gabor noise follows the general structure of that of sparse convolution noise [START_REF] Lewis | Algorithms for solid noise synthesis[END_REF][START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]]. We provide the annotated source code of the GPU implementation of our new noise in the supplemental material. In this section, we only highlight some important aspects of our implementation.

Coordinate Systems

To evaluate our filtered sliced solid Gabor noise (Eqn. 11), a point x, a unit normal n, and a filtering Gaussian G f are required. The point x and normal n are provided in (solid or 3D) texture space (tex), i.e., the space in which n (x; a, ω) is defined. The filtering Gaussian G f is provided to Eqn. 11 in tangent space (tan), i.e., the space determined by the 2D plane Π, but is defined in screen space (scr). The 2D plane Π is determined by xtex and ntex, which, together with two additional orthonormal vectors in Π, ttex and btex, determine a coordinate system for tangent space, and the mapping from texture space to tangent space Mtex→tan. The point xtex and normal ntex are obtained from their counterparts in object space (obj) and a (solid or 3D) texture placement matrix M obj→tex .

Continuity of the Local Coordinate System

We note that the local coordinate system formed by t, b and n does not have to be continuous over the surface. The tangents t and b do not affect the final filtered noise value, and therefore they can be chosen arbitrarily. We use an arbitrary normalized vector perpendicular to ntex for ttex, or the normalized version of ∂xtex/∂xscr,x, which lies in Π and is needed to compute Mscr→tex anyway.

Definition of the Filtering Gaussian

We define the 2D filtering Gaussian G f in a manner similar to [1989]. The filtering Gaussian is defined by its variance matrix in screen space, Σ f,scr , which typically is σ 2 f,scr I with σ f,scr = 1/2, and it is propagated to tangent space, in which G f should integrate to one, using the local affine approximation of the mapping from screen space to tangent space Mscr→tan. The filtering Gaussian in tangent space is thus G f,tan (xtan; 0, Σ f,tan ) with Σ f,tan = Mscr→tanΣ f,scr M T scr→tan . We obtain Mscr→tan on the GPU using the local affine approximation of the mapping of screen space to texture space, Mscr→tex, which is the matrix determined by the vectors ∂xtex/∂xscr,x and ∂xtex/∂xscr,y. These are easily obtained using built-in derivative functions (e.g., dFdx() and dFdy() in GLSL). Note that this definition also supports anisotropic screen-space filtering Gaussians (e.g., for motion blur).

Transforming the Phase-Augmented Gabor Kernel Slicing a 3D Gabor kernel (Eqn. 6) involves a transformation between texture space and tangent space, since the 3D kernel is defined in texture space but G f and the 2D sliced kernel are defined in tangent space. We note that transforming a phase-augmented Gabor kernel g (x; a, ω, φ) with a special orthogonal matrix M (i.e., a rotation) results in the phase-augmented Gabor kernel g x ′ ; a, ω t , φ with x ′ = M x and ω t = M ω.

Avoiding Matrix Inversions when Filtering Filtering a Gabor kernel (Eqn. 10) involves several matrix inversions. We reduce this number by exploiting the structure of the matrices (e.g., the covariance matrices Σ are symmetric and positive definite, and ΣG is a scaled identity matrix), and by simplifying the matrix expressions (see supplemental material for details).

Performance Our new random-phase solid Gabor noise runs at 43.74 and 28.83 megapixels per second for Fig. 4(a) and (b) (NVIDIA GeForce GTX 580 GPU,512×512,98,. We use a simple scene, since the performance figures are optimization-, implementation-and scenedependent. Despite the fact that the other noise functions have a reduced feature set, we also provide an indicative performance comparison. Our new noise is 1.8 -1.4× faster than Gabor surface noise (247 and 209 FPS), 2.4 -1.6× faster than wavelet noise (184 FPS), but up to an order of magnitude slower than Perlin noise (3050 FPS). Note that, although Perlin noise is faster than our noise, this extra performance cannot make up for the filtering quality: even with 256× stratified super-sampling, some of our "groundtruth" video's still exhibit aliasing (see e.g., supplemental material, Video fig1_c2.avi). Also note that the projection step in wavelet noise incurs a large performance penalty.

Results, Comparisons and Discussion

In Fig. 1, we compare Perlin noise [START_REF] Perlin | An image synthesizer[END_REF]], wavelet noise [START_REF] Cook | Wavelet noise[END_REF], surface Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]], and our new solid random-phase Gabor noise. We filter Perlin noise and wavelet noise using frequency clamping [START_REF] Norton | Clamping: A method of antialiasing textured surfaces by bandwidth limiting in object space[END_REF] (see supplemental material for details). We obtain ground-truth filtering using 256× stratified super-sampling. Our new noise function is the only one that supports high-quality filtering and does not introduce discontinuities at sharp edges. Please also refer to the videos in the supplemental material, which illustrate this more clearly. Note that we do not compare with [START_REF] Goldberg | Anisotropic noise[END_REF] since they do not support solid noise.

In Fig. 6, we compare a procedural solid texture generated with Perlin noise, wavelet noise, surface Gabor noise, and our new randomphase solid Gabor noise. We use a single scene (Fig. 6(a), a partial model of the Parthenon), a single texture (Fig. 6(b), a granite-like material), and two viewpoints (Fig. 6(a), red and green viewpoint). We compare the continuity at sharp edges (Fig. 6(c), view from the green viewpoint), and the filtering quality (Fig. 6(d), close-up from the view from the red viewpoint). We use the procedural texture model of Lagae et al. [2010b] (without the optional histogram matching), with parameters derived from exemplars and/or edited manually, and anisotropic noise octaves for the anisotropic examples. We disable lighting in order to avoid the introduction of lighting discontinuities. This comparison confirms the result of the previous one: Only our new noise function both supports high-quality filtering and does not introduce discontinuities at sharp edges. Please also refer to the images in the supplemental material, which illustrate this more clearly.

In Fig. 7, we show several solid procedural textures generated with our new random-phase solid Gabor noise. Our new noise function enables a variety of isotropic and anisotropic solid noise textures. Note the continuity across sharp edges (e.g., in Fig. 7(a), at the sharp edges of the chess pieces), and the high-quality anisotropic filtering (e.g., in Fig. 7(c), at the far end of the chess board). Perlin noise and wavelet noise do not support anisotropy. Surface Gabor noise does support anisotropy, but models a surface texture rather than a solid texture, which has a different appearance. Also note that, in contrast to surface Gabor noise, our noise does not require a continuous vector field over the surface. [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF] briefly mention solid Gabor noise as well. However, their solid noise does not support filtering; they simply evaluate 3D Gabor kernels at the surface.

It may seem that it would be possible to slice zero-phase rather than random-phase Gabor noise, and avoid the use of (random) phases altogether. However, this is not true, since slicing zero-phase Gabor noise introduces phase shifts in the sliced noise as well, and thus requires derivations similar to those in Sec. 3. Moreover, as pointed out in Sec. 3.4, the zero-phase-case is actually more complex than the random-phase case.

Conclusion

We have presented a new procedural noise function, random-phase Gabor noise, that supports continuity across sharp edges, highquality anisotropic filtering, and anisotropy. We have achieved this by augmenting Gabor noise with phase, which allowed us to explicitly model the slicing of the 3D Gabor kernels, which in turn allowed us to filter the resulting 2D Gabor kernels. We believe that Gabor noise should always be used with random phases, since it both strengthens and simplifies the model.

Interesting opportunities for future work include exploring volumetric filtering of solid Gabor noise (e.g., for hyper-texture, when ray marching through solid noise), further exploring anisotropy in the context of solid noise, and designing user interfaces for interacting with anisotropic solid noise.

Figure 1 :

 1 Figure 1: Existing noise functions either introduce discontinuities of the solid noise at sharp edges, which is the case for wavelet noise (b.1) and Gabor noise (c.1), or result in detail loss when anti-aliased, which is the case for Perlin noise (a.2) and wavelet noise (b.2). We present a new noise function that preserves continuity over sharp edges (d.1) and supports high-quality anti-aliasing (d.2).

Figure 2 :

 2 Figure 2: Slicing and projecting a band-pass noise. Slicing in the spatial domain (solid green line) (a) corresponds to projecting in the frequency domain (integration in the direction of the dashed green line, resulting in the red curve) (c), which does not result in a band-pass power spectrum. Projecting in the spatial domain (integration in the direction of the dashed green line) (b) corresponds to slicing in the frequency domain (solid green line) (d), and introduces discontinuities at sharp edges.

Figure 3 :

 3 Figure 3: Slicing the Gabor and the phase-augmented Gabor kernel: (a) Slicing the 2D Gabor kernel; (b) Slicing the 2D Gaussian; (c) Slicing the 2D harmonic; (d) Slicing the 2D phaseaugmented Gabor kernel. Note the phase shift φ s introduced by slicing.

Figure 4 :

 4 Figure 4: Different kinds of solid random-phase Gabor noise: (a) anisotropic; (b) isotropic; (c) hybrid anisotropic/isotropic.In its most basic form (Eqn. 7), solid random-phase Gabor noise is an anisotropic solid noise parameterized by a (3D) frequency ω and a bandwidth a (Fig.4(a)). Other kinds of solid noise are obtained depending on how the parameters {ai} and {ω i } vary for different kernels. We obtain an isotropic solid noise parameterized by a (1D) radial frequency f and a bandwidth a (Fig.4(b)) using ωi,r = f , ω i,θ ∼ U [0, 2π] and ω i,φ ∼ cos -1 {U [-1, 1]}. We obtain an interesting hybrid anisotropic/isotropic solid noise (Fig.4(c)) similarly but with ω i,φ = 0. A greater variety of solid noise patterns could be obtained by compositing solid noise patterns using widgets, in[START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]], or by using spatially varying solid noise parameters, as in[START_REF] Lagae | Improving Gabor noise[END_REF]].

Figure 5 :

 5 Figure 5: The Slice-Projection Theorem for Stochastic Processes: The power spectrum of X is SXX ; therefore, the power spectrum of Y is SY Y .phases. Intuitively, this is because the random phases ensure invariance w.r.t. slicing: since the phases of the Gabor kernels are random, the phase-shifted phases of the sliced kernels are random as well. In other words, the random phases allow us to both strengthen and simplify the Gabor noise of[START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF].

Figure 6 :

 6 Figure 6: Comparison of a procedural solid texture generated with Perlin noise, wavelet noise, surface Gabor noise and our new randomphase solid Gabor noise, using the procedural texture model of Lagae et al. [2010b]: (a) Scene, partial model of the Parthenon; (b) Texture, granite-like material; (c) Continuity across sharp edges, green viewpoint (no discontinuity is better); (d) Filtering quality, red viewpoint (more detail towards the right is better). Our new noise function is the only one that supports high-quality filtering and does not introduce discontinuities at sharp edges.

Figure 7 :

 7 Figure 7: Solid procedural textures generated with our new random-phase solid Gabor noise, using the procedural texture model of Lagae et al. [2010b]: (a) Granite chess pieces; (b) Granite dancer statuettes; (c) Wooden chess board. Our new noise function enables a variety of isotropic and anisotropic solid noise textures. NVIDIA (Academic Partnership Program). Ares Lagae is a Postdoctoral Fellow of the Research Foundation -Flanders (FWO), and acknowledges K.U.Leuven CREA funding (CREA/08/017).

Acknowledgements

We would like to thank the anonymous reviewers; Carles Bosch, Fredo Durand, Bruno Galerne, Ian Jermyn, Sylvain Lefebvre, Javier Taibo (maya2osg) and Celine Vens; and Autodesk (software donation) and


Julie Bertrand
email: julie.bertrand@inserm.fr

Céline M Laffont

France Mentré

Marylore Chenel

Emmanuelle Comets

Development of a complex parent-metabolite joint population pharmacokinetic model

Keywords: Nonlinear mixed effects modeling, First-pass metabolism, Back-transformation mechanism, Estimation algorithms, Genetic covariate

This study aimed to develop a joint population pharmacokinetic model for an antipsychotic agent in development (S33138) and its active metabolite (S35424) produced by reversible metabolism.

Because such a model leads to identifiability problems and numerical difficulties, the model building was performed using the FOCE-I and the SAEM estimation algorithms in NONMEM and MONOLIX, respectively. Four different structural models were compared based on Bayesian information criteria. Models were first written as ordinary differential equations systems and then in closed form (CF) to facilitate further analyses. The impact of polymorphisms on genes coding for the CYP2C19 and CYP2D6 enzymes respectively involved in the parent drug and the metabolite elimination were investigated using permutation Wald test. The parent drug and metabolite plasma concentrations of 101 patients were analyzed on two occasions after 4 and 8 weeks of treatment at 1, 3, 6 and 24 hours following daily oral administration. All configurations led to a two compartment model with back-transformation of the metabolite into the parent drug and a first-pass effect. The elimination clearance of the metabolite through other processes than back-transformation was decreased by 35 [9-53] % in CYP2D6 poor metabolizer.

Permutation tests were performed to ensure the robustness of the analysis, using SAEM and CF.

In conclusion, we developed a complex joint pharmacokinetic model adequately predicting the impact of CYP2D6 polymorphisms on the parent drug and its metabolite concentrations through the back-transformation mechanism.

INTRODUCTION

For most drugs undergoing biotransformation, the process results in the formation of a more polar metabolite that is pharmacologically inactive and is eliminated more rapidly than the parent drug [START_REF] Shargel | Applied Biopharmaceutics & Pharmacokinetics[END_REF]. For some drugs however, the metabolite may be pharmacologically active and/or produce adverse effects. Joint pharmacokinetic (PK) modelling accounts for the uncertainties in the data and allows feedback from the metabolite data to the parent drug data to influence the estimation. Thus, more and more studies now use this approach [START_REF] Younis | Enterohepatic recirculation model of irinotecan (CPT-11) and metabolite pharmacokinetics in patients with glioma[END_REF][START_REF] Panhard | Population pharmacokinetic analysis for nelfinavir and its metabolite M8 in virologically controlled HIVinfected patients on HAART[END_REF][START_REF] Feng | Population pharmacokinetic analysis for risperidone using highly sparse sampling measurements from the CATIE study[END_REF]. Also, joint modelling allows to correctly evaluate and predict the impact of drug-drug interactions, and in the case of reversible metabolic systems, it is the only way to properly assess covariate effects.

However, such models can rapidly gain in complexity and present parameter identifiability problems and numerical difficulties in terms of estimation. To solve the issue of structural identifiability, it is important to identify the parameters or "apparent" parameters that can be estimated [START_REF] Panhard | Population pharmacokinetic analysis for nelfinavir and its metabolite M8 in virologically controlled HIVinfected patients on HAART[END_REF], and for the sake of parameter interpretation or covariate analysis, it might be needed to make some assumptions on the parameters (e.g. fixing one parameter to a given value) [START_REF] Feng | Population pharmacokinetic analysis for risperidone using highly sparse sampling measurements from the CATIE study[END_REF]. Also, instability during the estimation is likely to result in solutions at local minima which may lead to biased parameter estimates and potentially wrong conclusions. Furthermore, numerical issues are likely to arise depending on the algorithms used for estimation in nonlinear mixed effect models (NLMEM). The First Order Conditional Estimation with Interaction (FOCE-I) algorithm [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF] which is implemented in the NONMEM software [START_REF] Beal | NONMEM User's Guides[END_REF], has been shown to encounter numerical difficulties even on simple single response models [START_REF] Girard | A comparison of estimation methods in nonlinear mixed effects models using a blind analysis[END_REF][START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF]. Yet, more robust alternatives have been proposed that avoid simplifying the equation for the likelihood, such as the Stochastic Approximation Expectation Maximization (SAEM) algorithm [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] implemented in the MONOLIX [START_REF] Lavielle | MOdèles NOn LInéaires à effets miXtes)[END_REF] and S-ADAPT [START_REF] Bauer | S-ADAPT/MCPEM User's Guide Version 1.56[END_REF] software. The NONMEM software remains though the most popular tool in population PK analysis because of its superior flexibility in the pharmaceutical field with tools such as NMTRAN. MONOLIX was only recently provided with a similar model translator called MLXTRAN. Of note, the SAEM algorithm was recently implemented in NONMEM software version 7.2.

During its development, an innovative antipsychotic agent from SERVIER research, the S33138 [START_REF] Millan | 3aS,9bR)-8-cyano-1,3a[END_REF], was shown to produce a metabolite which was also an active compound, called S35424, through hydrolytic cleavage by hydrolases, as assessed in human microsomes. Preclinical studies in rat and monkey supported the existence of a back-transformation of the metabolite into the parent drug, and this was further confirmed in a microdose study in man. Such reversible metabolism explained why both compounds showed similar terminal plasma half-life although the parent compound only slightly accumulates while the metabolite accumulates to a much larger extent (AUC accumulation ratio of 3). In addition to hydrolases, other metabolism pathways were identified for the parent drug involving CYP3A4 and to a lesser extent CYP2C19. The active metabolite was mainly metabolized through CYP2D6. Both CYP2C19 and CYP2D6 are encoded by highly polymorphic genes and these polymorphisms are known to have an impact on the course of many therapeutic drugs (13). Therefore, the effects of CYP2D6 and CYP2C19 polymorphisms on the PK of the parent compound and the metabolite were investigated in a phase II study using NLMEM. Drug plasma concentrations profiles were documented in patients at two occasions for the parent drug and its metabolite along with genotypes for five polymorphisms of the CYP2D6 gene and two polymorphisms of the CYP2C19 gene. In previous works, we have shown through simulations that asymptotic tests to detect a gene effect in NLMEM require a correction for type I error inflation on designs with unbalanced genotypes and/or including a small number of subjects [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF][START_REF] Bertrand | Pharmacogenetics and population pharmacokinetics: impact of the design on three tests using the SAEM algorithm[END_REF]. This slight inflation can be handled by permutation tests [START_REF] Manly | Bootstrap and Monte Carlo Methods in Biology[END_REF]. However such computing intensive approaches require a fast estimation method, all the more when the structural and variability models become complex.

The aim of the present work was thus to develop a joint population PK model for the antipsychotic and its metabolite after oral administration based on the data of the phase II study and to test for genetic effects. Because of the model complexity (reversible metabolism), the building of the model was performed using FOCE-I in NONMEM and SAEM in MONOLIX, in parallel. Furthermore, in order to alleviate the computational burden of permutation tests in the analysis of the genetic effect of CYP2D6 and CYP2C19 polymorphisms, models were encoded in ordinary differential equations (ODE) system and closed form solutions (CF). Finally, both internal and external model evaluations were performed, the latter using the data from a phase I study. It is noteworthy that external evaluation provides the most stringent method for assessing the predictive ability of the developed model [START_REF] Brendel | Are population pharmacokinetic and/or pharmacodynamic models adequately evaluated? A survey of the literature from 2002 to 2004[END_REF].

MATERIALS AND METHODS

Pharmacokinetic studies

The data used for model building came from a pilot, phase II, international, multicentre, randomized, double-blind, parallel-group, study assessing the effect of the above-mentioned antipsychotic drug compared to a gold-standard medication. In this study, one hundred and twenty patients were randomly allocated to 4 groups receiving for 8 weeks either 5, 10 or 20 mg of the novel antipsychotic or the gold standard medication (risperidone). Both treatments were administered orally and once a day in the morning before breakfast. Four blood samples were collected on two occasions, i.e. four and eight weeks after treatment initiation (W4 and W8). The sampling times were set empirically as the following: prior to drug administration and then 1, 3 and 6 h after drug administration. On the day of the sampling, the date and exact time of each blood sample were reported along with the date and exact time of the drug administration as well as those of the previous administration. Those exact times were used for the modelling.

For the external model evaluation, we used the data from a phase I randomized, double-blind tolerance study versus placebo with repeated increasing oral doses. Thirty healthy male volunteers were randomly divided into three groups: in each group, eight subjects received repeated once-a-day oral administrations of the parent drug at the dose of 10, 20 or 30 mg respectively, while two other subjects received the placebo. In this study, blood concentrations were documented at steady-state after 14 days of treatment with the following sampling times: prior to drug administration and then 0.33, 0.66, 1, 1.5, 2, 3, 4, 8, and 12 hours after drug administration. The exact times of blood collection were also recorded.

Concentration measurements and genetic polymorphisms

In both studies, the plasma concentrations of the antipsychotic and its active metabolite were determined using a validated method involving solid phase extraction followed by reverse phase Liquid Chromatography with Mass Spectrometry -Mass Spectrometry detection (LC/MS-MS).

The intra/inter-day accuracy were 97.8/101.7% for the parent drug and 93/101.3% for the metabolite, the intra/inter-day precision were 4.7/5.6% for the parent drug and 6.1/8.6% for the metabolite and the recovery was 71% for the parent and 65.4% for the metabolite. Typical retention times of the chromatography were about 1.6 min for S33138 and the internal standard and 2.0 min for S35424, respectively. The limits of quantification were 0.56 and 0.44 nmol.L -1 for the parent drug and its metabolite, respectively.

In the phase II study, blood samples were taken at the selection visit in order to determine the patient's genotype of CYP 2D6*3, *4, *6, *7, and *8 alleles and of CYP 2C19*2 and *3 alleles while in the phase I study, a phenotyping test (dextrometorphan administration) was performed at the selection visit to ensure that all the persons involved were CYP2D6 intermediate or extensive metabolizers. All samples were stored at -20 • C and blood samples of included patients were sent on a regular basis to a central laboratory for analysis.

Joint pharmacokinetic model

Here, we present the analysis of concentration-time profiles from patients treated with the novel antipsychotic only (excluding patients treated with the reference medication). Concentrations were expressed in nmol.L -1 for the joint modelling of the parent drug and its metabolite; molecular weights were respectively: 319.4 g.mol -1 and 361.4 g.mol -1 . Data below the limit of quantification as well as concentrations considered as non reliable given time or dosing information were treated as missing data and were excluded from the analysis.

Four different structural models of increasing complexity were investigated: their structural representations are shown in Figure 2 along with the definition of related parameters. The first model, on the left end of Figure 2, is a two compartment model where the dose is absorbed into the parent compartment at a rate K ap and the parent drug (S33138) is either eliminated from the system with a clearance CL po or is transformed with a clearance CL pm into a metabolite (S35424) which is eliminated from the system with a clearance CL mo . Panhard et al. used such a model to jointly analyze nelfinavir and M8 concentrations [START_REF] Panhard | Population pharmacokinetic analysis for nelfinavir and its metabolite M8 in virologically controlled HIVinfected patients on HAART[END_REF]. In the second model from the left on Figure 2, S35424 can in addition be transformed back to the parent with a clearance CL mp . The existence of a back-transformation of S35424 into S33138 was supported by preclinical studies in rat and monkey and further confirmed in a microdose study in man. This microdose study included five healthy volunteers in a cross-over design where they received an oral dose of S33138 plus an additional intravenous dose of C 14 S35424 or C 14 S33138 depending on the period with two groups of sequence composed of three and two subjects, respectively. Substantial concentrations of C 14 S33138 were measured after administration of C 14 S35424. Moreover, back-transformation of the metabolite into the parent drug is a well known process for numerous amines. A similar model was used to describe the increased hydrolysis of tesaglitazar metabolite into tesaglitazar via biliary circulation in renally impaired subjects [START_REF] Hamrén | Mechanistic modelling of tesaglitazar pharmacokinetic data in subjects with various degrees of renal function-evidence of interconversion[END_REF]. In the third model from the left on Figure 2, a first-pass effect has been considered: the dose not only enters into the parent compartment at a rate K ap but also enters into the metabolite compartment at a rate K am . This model, although less physiological, is parsimonious because the fraction of dose transformed into the metabolite is driven by the two absorption rate constants. The last model, on the right end of Figure 2, includes a dose apportionment independent of the rate constant values (K ap and K am ) with a fraction F p of the dose leading to the parent and a fraction 1-F p leading to the metabolite prior to reach the plasma. Such presystemic formation of an interconversion metabolite is mentioned in the extensive review on reversible metabolic systems from Cheng and Jusko [START_REF] Cheng | Pharmacokinetics of reversible metabolic systems[END_REF]. For identifiability purposes, the fraction of dose available after absorption (f ) was set to 1.

Also, for all these models, it is noteworthy that all PK parameters cannot be estimated in the present study, since only an oral administration of the parent drug was performed. Estimating all PK parameters requires that both the parent drug and the metabolite be given by the intravenous route in addition to the oral administration of the parent compound [START_REF] Cheng | Pharmacokinetics of reversible metabolic systems[END_REF]. It results that only "apparent" parameters can be estimated. The following parameters were identified for the last model:

V p /F p , V m /(1-F p ), E 1 = k po +k pm , E 2 = k mo +k mp , k pm F p /(1-F p ) and k mp (1-F p )/F p . An
alternative parameterization is given in Appendix. As our objective was here to assess the impact of genetic covariates on specific clearances (CL mo for CYP2D6 and CL po for CYP2C19 ), it was necessary to make an assumption on the parameters to estimate all clearances separately (and not global clearances divided by the corresponding volume). Thus, we chose to set V m equal to V p . In this respect, F p , V p and the clearance estimates presented here are tagged by a * as a reminder of their reliance on the assumption made on volumes. All drug transfers between compartments were modelled as linear processes.

Population parameter estimation was performed in parallel with the FOCE-I and SAEM algorithms implemented in the NONMEM and MONOLIX softwares, respectively. The first algorithm performs an approximation of the model with a first-order linearisation around the individual predictions of the random effects. In contrast, the second algorithm is a stochastic version of the well-known Expectation-Maximization algorithm where the individual random effects are the missing variables. In the expectation step, the individual parameters are simulated using a Monte-Carlo Markov Chain approach and then used to compute a stochastic approximation of the conditional expectation of the complete log-likelihood which is the log-likelihood of the complete data, i.e the observations and the imputed individual parameter estimates at the current iteration of the SAEM algorithm. Then , the complete log-likelihood is maximized to obtain the updated estimates of the population parameters. In NONMEM, the convergence of the FOCE-I algorithm is achieved when all parameters have the required number of significant digits, while in Expectation Maximization-like methods the definition of convergence is usually left to the user discretion.

In MONOLIX additionally to the visual inspection of convergence graphics, stopping rules are available which are based on the absence of decrease of the complete log-likelihood sequence during the stochastic step and small variability between subsequent population parameter estimates and estimates of the complete log-likelihood during the cooling step.

In both NONMEM and MONOLIX, the Fisher Information Matrix and the loglikelihood of the model were obtained by linearization of the model around the predictions of the individual parameters. In NONMEM, we used the OPTION UNCONDITIONAL in the covariance step to overcome the convergence issue as well as the OPTION MATRIX=R and MATRIX=S when needed.

The models were encoded through ODE as well as using the corresponding CF solution derived using the Laplace transform approach (19, 20) (see appendix). For the models encoded with ODE system, data were fitted using the ADVAN5 routine in the NONMEM software version 7.2 and MLXTRAN with the STIFF option in MONOLIX version 2.4. For the models with a dose apportionment, we used the ADVAN6 routine because of an issue in the ADVAN5 routine to be fixed in the next version of NONMEM. For the models encoded in CF solution, data were fitted using the PRED routine in the NONMEM software version 7.2 and the model building function of MONOLIX version 2.4.

In NONMEM, the number of digits required for convergence (SIG) was set to its default value of 3 for ADVAN5 and to 3 or 2 for ADVAN6 with a fortran 95 compiler. With ADVAN6, we set the number of significant digits for the predicted values (TOL) to 6 and the number of significant digits for the objective function (SIGL) to 6 in the estimation step and TOL=6 with SIGL=6 in the covariance step, following guidance from [START_REF] Beal | NONMEM User's Guides[END_REF]; with ADVAN5 the TOL option is not required.

In MONOLIX, the algorithm settings were left to the default values; the maximal numbers of stochastic (K1) and cooling (K2) iterations were set to 500 and 200 respectively, with use of automatic stopping rules and only one markov chain (nmc).

The structural model was determined on the data collected at W4 only and all the systems were considered at steady state with a 24 hours dosing interval. We used an exponential model for the between subject variability to ensure positivity of the individual PK parameters and a logit model to force 0 ≤ F p ≤ 1. Along the structural model selection, random effects for all parameters but f were assumed to follow a normal distribution with a diagonal variance matrix.

A combined error model with an additive component a and multiplicative coefficient b was used for the parent drug and the metabolite. Model selection was performed using the Bayesian information criteria ( 21): BIC = -2L + P pop log(N ), where L is the loglikelihood of the model, P pop is the number of population parameters which includes the fixed effects, and the variance components and N the number of subjects. BIC allows to compare both nested and non nested models such as those including F p or not. In order to compare the NONMEM and MONOLIX output, we retrieved -2L for each model by adding n tot × log(2π) to the objective function estimate, n tot being the total number of observations. Also in the present work we used N rather than n tot in the BIC penalizing term, as the observations are dependent within each subject [START_REF] Faes | The Effective Sample Size and an Alternative Small-Sample Degrees-of-Freedom Method[END_REF][START_REF] Visser | Effective Sample Size and the Bayes Factor[END_REF].

Once the structural model was selected, parsimonious error models (additive and proportional) were tested using the likelihood ratio test (LRT). To capture clearance and volume correlations due to bioavailability variation, between-subject variance was estimated on f . Next, between-subject variance (ω 2 ) nullity was tested using a LRT with, as a reference distribution, a mixture of a χ 2 distribution with 0 degree of freedom and χ 2 distribution with 1 degree of freedom [START_REF] Stram | Variance components testing in the longitudinal mixed effects model[END_REF]. Finally, data at W8 were added to the data set and forward inclusion of within-subject variances (γ 2 ) was performed on parameters with non-null between-subject variance.

Covariate model

A linear dose effect was systematically investigated on f and F p parameters using the Wald test.

The dose was analyzed as a continuous covariate with 10 mg as the reference dose, so that f and F p of subject i at occasion k were predicted as described in equations ( 1) and [START_REF] Younis | Enterohepatic recirculation model of irinotecan (CPT-11) and metabolite pharmacokinetics in patients with glioma[END_REF].

f ik = f e β f,D (D i -10)
(1)

F p,ik = Fp 1-F p e β Fp,D (D i -10) 1 + Fp 1-F p e β Fp,D (D i -10) (2)
Also, linear effects of the CYP2D6 and CYP2C19 polymorphisms were investigated following a forward selection based on the Wald test with a significance threshold at 5%. Both genetic covariates were analyzed by means of a phenotypic binary categorization; poor metabolizer (PM)

versus the reference class the extensive metabolizer (EM), so that the parameter θ of subject i at occasion k was predicted as described in equation ( 3) where β θ,CY P * was non-null for patients CYP* EM and 0 otherwise.

θ ik = θ e β θ,CY P * (3)
The classification was performed as follows: carriers of two rare allele were classified as PM [START_REF] Chou | Comparison of two CYP2D6 genotyping methods and assessment of genotype-phenotype relationships[END_REF].

In order to correct the inflation of the Wald type I error shown to occur on unevenly distributed genotypes ( 14), final p-values were assessed using permutations [START_REF] Manly | Bootstrap and Monte Carlo Methods in Biology[END_REF]. The null hypothesis of the permutation test is that the mean value of the parameter θ would be the same whatever the CYP2D6 and CYP2C19 status while the alternative is that the mean value would be different between the poor and extensive metabolizers. More specifically, 1000 data sets were generated by permuting the rows of the covariate matrix from the original data set. For each covariate, one Wald statistic, W obs , was estimated from the original data and one Wald statistic, W perm , was estimated from each of the 1000 data sets. Thus, we obtained j=1,...,1000 W perm j . The permutation p-value was the proportion : (card(W perm j ≥ W obs )+1)/(1000+1). As the patients were sampled on two occasions, we also computed a metric representing the genetic component 26), for each model parameter with non-null withinsubject variance. This component gets closer to 100% when the between-subject variance for the parameter under study is larger than its within-subject variance, so that the variability for this parameter is more likely to be explained by a genetic covariate.

of variability R GC = max 0, 1 -γ 2 ω 2 × 100 (
As mentioned previously, the covariate analysis was performed with the assumption V m = V p .

As the parameter estimates are dependent on this assumption, we have also derived apparent parameters that are independent of any assumption made:

V p /F p , V m /(1 -F p ), CL ptot /F p ,
CL mtot /(1 -F p ), CL pm /V m and CL mp /V p (see appendix). The fixed effects and between-subject standard deviation of these apparent parameters were derived by simulation using our final parameter estimates.

Model evaluation

An internal model evaluation was performed on both occasions as well as an external model evaluation which only considered one occasion.

For the internal model evaluation before the addition of the genetic effect, we performed visual predictive check plots where the 90% confidence intervals around the 5 th , 50 th and 95 th prediction percentiles from 250 simulated data sets were overlaid to the 5 th , 50 th and 95 th percentiles of the observed data binned using the theoretical sampling times [START_REF] Karlsson | 17th Population Approach Group in Europe[END_REF]. Then after the addition of the genetic effect, we computed the normalized prediction distribution errors (npde), i.e. the observation percentiles within the empirical distribution obtained from the model simulations, decorrelated and normalized using the inverse function of the normal cumulative density function [START_REF] Brendel | Evaluation of different tests based on observations for external model evaluation of population analyses[END_REF].

For the external evaluation of the model after the addition of the genetic effect, we computed the npde and performed a visual predictive check plots where the 5 th , 50 th and 95 th predicted percentiles from 250 simulated data sets were overlaid to the observed 5 th , 50 th and 95 th percentiles of the observed data binned using the theoretical sampling times [START_REF] Karlsson | 17th Population Approach Group in Europe[END_REF].

For the internal and the external model evaluation based on npde, 1000 data sets were simulated and we used the R npde package (29) for the calculations.

RESULTS

Data

The and 53.5% were men. Among the 101 patients: 12 patients were classified as CYP2D6 PM and 2 were classified as CYP2C19 PM. Four of the CYP2D6 PM patients received a dose of 5 mg, 6 received a dose of 10 mg and 2 received a dose of 20 mg. The two CYP2C19 PM patients received a dose of 5 mg and 20 mg, respectively. Genotype information was not available in 2 patients.

On W4, 355 and 358 concentrations were measured for the parent drug and the metabolite respectively in 97 patients, 81 of whom had a complete profile of four samples for both compounds.

On W8, 271 and 268 concentrations were measured for the parent drug and the metabolite respectively in 71 patients, 51 of whom had a complete profile for both compounds. Sixty-seven patients had concentration-time profiles at both W4 and W8. Only 8 concentrations out of 1252 (less than 1%) were below the limit of quantification and were discarded from the analysis.

Figure 2 displays the observed concentrations for the parent drug and the metabolite on both normal and semi-log scales versus the time at W4 and W8, revealing a much larger accumulation for the metabolite than for the parent drug.

In the external evaluation data set, the average age was 28 years , the average weight 72.5 kg [63.0-99.3] and all patients were men. Eight patients had a dose of 10 mg, 8 a dose of 20 mg and 7 a dose of 30 mg. After two weeks of repeated oral doses, 203 concentrations were measured for the parent drug and the metabolite in 23 healthy volunteers, 20 of whom had a complete profile of ten samples for both compounds. No concentrations were below the limit of quantification.

Structural and variability models

Using FOCE-I in NONMEM and SAEM in MONOLIX, convergence was achieved for all models with both coding conditions (ODE or CF). However with FOCE-I in NONMEM, standard error estimates could only be obtained for the models with dose apportionment in ODE.

Table I reports the BIC along with the error model variance component estimates as well as

the computation time for the four structural models investigated plus the model with K ap = K am , using both software and codings on data at W4. Confirming pre-existent knowledge, the addition of the back-transformation mechanism greatly improved the model with a 100 units drop in BIC. The addition of a first-pass effect led to a 40 units decrease in BIC further lowered with the inclusion of an additional F p term for the dose apportionment. Estimates of error model components were quite close across algorithms and coding with the exception of the additive coefficient for the parent drug. Yet, the relative standard error (RSE) obtained with FOCE-I in NONMEM and/or MONOLIX for this parameter was always very large (data not shown). For the most complex model, computation times dramatically decreased when using CF instead of ODE irrespective of the software (53.8 h to 22.2 min using NONMEM and 2.42 h to 0.6 min using MONOLIX).

The fixed effect estimates were similar when using ODE or CF for the same model and estimation algorithm, as shown in Table II. The median and maximal relative difference on clearances were respectively 0 and 47% with FOCE-I and 3 and 45% with SAEM in MONOLIX. A proportional and a combined error model were selected for the parent drug and its metabolite, respectively, with estimates for the multiplicative component about 30% and 6%. Of note in the model including the f parameter, the BIC was reduced by 100 units using a diagonal-covariance matrix for the random effects compared to a full-matrix where the smallest and largest correlation estimates were -0.09 and 0.71, respectively. Then, estimation of between-subject variance was found to improve BIC for all parameters except the fraction of dose that escaped first-pass effect (F p *), CL mp * and CL pm *. Once included the data at W8, the absorption constant rate had to be fixed to its estimate on data at W4 for stability purposes. Only the within-subject variance of the CL po * was significantly different from 0.

Covariate model

It appears that the f * and F p * were 10 and 22% higher for a 5 mg dose and 19 and 33% lower for a 20 mg dose respectively (both p-values< 10 -3 ). The population parameter estimates of the basic model not including the genetic covariate and their relative standard error (RSE) are given in Table IV, for the 101 patients of the phase II study. Figure 4 represents the visual predictive check (VPC) plots obtained for each compound at both occasions and for the three doses. The predictions from the model described adequately the observed high and median concentration profiles of both molecules for the three doses. However, the model seems to predict less well the low concentration-time profiles, with a systematic overprediction at time 1h for the parent together with an underprediction of the ensuing times for the 5 mg dose on both occasions and the 10 mg dose at W8. For the metabolite, the VPC plots indicate a misfit for the low concentration-time profiles for the 10 mg dose at W4. It is noteworthy though that some of these misfits might be explained by some atypical individual profiles (see figures 3c and 3d), which drive the lower observed percentile estimates at doses 5 and 10 mg.

The distribution. Yet at time 1h on both occasions, the npde showed a deviation which was previously observed on the visual predictive check plots and the corresponding global tests p-values given in the output of the npde calculation were below 0.001. The bias at time 1h might be explained by the existence of a few patients with a slower absorption whose S33138 profile is not accurately described by the model. For the metabolite, the npde appeared randomly distributed around 0 and within the required boundaries but the corresponding global tests p-values were also below 0.001. Figure 7 and 8 present respectively the npde and a visual predictive check plots from the model including the genetic covariate versus time for the parent drug and the metabolite in the external evaluation data set. The npde for the parent drug tended to be below zero, especially in the first times but remained above the 5 th percentiles of the normal distribution. Also, the colored area formed by the 90% interval was somewhat narrower than expected and the corresponding global tests p-values were below 0.001. Similarly on the visual predictive check plots, the 5 th and 95 th predicted percentiles were respectively lower and higher compared to the observed percentiles, especially for the 95 th percentiles of the parent drug at the 20 mg dose. Nevertheless the observed and predicted median were satisfactorily close.

In conclusion, the final model provided sensible predictions of the observations despite some difficulties in capturing the lower concentration profiles in patients and an overprediction of the variability in healthy volunteers. Due to these misspecifications, the model was rejected based on the npde global test. However, the heterogeneity in absorption at 1h could be tackled in future analyses using mixture models to account for the diet or other unknown factors and healthy volunteers are known to be a less heterogeneous population than patients. This model, involving back-transformation and first-pass effect, had though to be parsimonious given the study design.

DISCUSSION

In the present work, we compared four different structural models to describe the concentrationtime profiles of a novel antipsychotic and its active metabolite obtained at two occasions. We used the FOCE-I algorithm in NONMEM and the SAEM algorithm in MONOLIX to fit the data using models encoded both in ODE and CF solutions. With the final selected model, we investigated the effect of the metabolizer status as defined by CYP2D6 and CYP2C19 polymorphisms on the PK of the parent drug and its metabolite.

The PK of both compounds was reasonably well described by a model including a back-transformation mechanism and a first-pass effect with an additional term for the dose apportionment. Cheng and Jusko ( 18) defined the acetylation of the metabolite to form back its acetylated parent as a known process for numerous amines, in their extensive review on reversible metabolic systems. They also mentioned the occurrence of presystemic formation of an interconversion metabolite and the consequences in determining the rates and extents of drug absorption and presystemic formation of an interconversion metabolite. In addition, this reference work lists several formulas to derive PK parameters for compounds with interconversion metabolite(s), based on AUC ratio and/or numerical deconvolution after a separate analysis of the parent and/or metabolite profiles. Here, we chose to model both compounds jointly, to provide PK parameter fixed effects and variance estimates in patients. In the present study, the backtransformation process provided a mechanism-based explanation for the long terminal plasma half-life observed for the parent drug with reduced accumulation at steady-state. Actually, the terminal plasma half-lives derived from our model for the parent drug and metabolite were similar, yet the parent drug showed an effective half-life of 8.5 h i.e. 4 times faster than that of the metabolite (32.7 h) [START_REF] Rowland | Clinical pharmacokinetics. Concepts and Applications[END_REF]. The first-pass effect, on the other hand, allowed to capture the early bump observed in the metabolite data. Auclair et al. also showed that adding a first-pass effect in their model allowed to fit metabolite concentrations that appeared quicker or at the same time as the parent drug [START_REF] Chou | Comparison of two CYP2D6 genotyping methods and assessment of genotype-phenotype relationships[END_REF]. Such a first-pass effect through amidases in the gut and the liver has already been observed in more dramatic proportions for the experimental anti-convulsant related to lidocaïne, D2624 [START_REF] Martin | Pharmacokinetics and metabolism of the novel anticonvulsant agent N-(2,6-dimethylphenyl)-5-methyl-3isoxazolecarboxamide (D2624) in rats and humans[END_REF]. In our final model, the rate of appearance of the metabolite, K am was set to the absorption rate of the parent drug K ap . We explain this finding by the fact that the absorption rate of the parent drug would be rate limiting due to the very quick formation of the metabolite. Duffull et al. also reported this phenomenon for ivabradine and its metabolite in [START_REF] Duffull | A pharmacokinetic simulation model for ivabradine in healthy volunteers[END_REF].

Thorough mathematical investigations were required to identify the apparent PK parameters which are presented here for the final and more complex model (first-pass effect with backtransformation of the metabolite) in the appendix of the article. Yet, interpretation of apparent PK parameters is not always straightforward, and they may not allow proper covariate analysis. For instance, the effect of CYP2D6 was expected on CL mo only, given prior knowledge, with no impact on CL mp so using the apparent parameter CL mtot /(1 -F p ) would not have been appropriate. For this reason, we have decided to make an assumption on volumes in order to be able to estimate all clearances separately. By convenience, we chose to set the volumes equal as this allowed estimation of between and within-subject variability on parent and metabolite volumes. This assumption was judged reasonable since the ratio of 1 is in the range of volumes ratios estimated from the microdose study and gives estimates for the bioavailability and the percentage of dose undergoing back-transformation close to those obtained from the microdose study where clearances of all processes were identifiable. The microdose study gave a median [range] volume of distribution of 78 L [60-287] for the parent and 48 L for the metabolite, with a volume ratio of 2.4 [0.5-6].

It is obvious that Fp, Vp and clearance estimates depend on the assumption we made on the ratio of volumes but setting the volumes ratio to 2.4 in the final model did not markedly impact the results of the analysis apart a significant higher estimate for the percentage of dose undergoing back-transformation (10 % instead of 5 %) which is less in agreement with the microdose study results.

In silico measures of log P=2.29 and 2.36, log D (pH7.4)=1.24 and 1.32 and pKa1=8.46 and 8.46 (tertiary amine) between respectively the parent and the metabolite were similar with the exception of pKa2=12.43 (amide) and 4.73 (primary amine). This information neither supported nor infirmed the assumption on the ratio of volumes. Moreover, the volume of distribution also depends on the compound affinity for transporters or tissues as well as its binding ability to plasma proteins. In the present analysis, models with saturable elimination were not investigated as neither the prior knowledge on the compounds nor the goodness of fit plots suggested the existence of such a process. Also, models with a central and peripheral compartment for the parent and metabolite were not evaluated because of the sparse study design which would have not enabled robust estimation.

Here, two estimation algorithms implemented in two softwares were used for data analysis.

As regulatory authorities encourage the use of new estimation algorithms, this work addresses a current and relevant issue with a comparison of both software on such a complex model and using both codings. The SAEM algorithm has been recently implemented in NONMEM version

7. Yet, we chose its implementation in MONOLIX because it has been thoroughly evaluated compared to the more recent NONMEM one. Both algorithms led to the selection of the same structural model whatever the coding approach. However, substantial differences between the standard deviation estimates were observed across codings for the same model and same estimation algorithm. An extensive comparison of both softwares performances has been already performed and is currently available on the MONOLIX website (http://software.monolix.org/sdoms/ software/index.php?/evaluation.html). However, the models investigated were less advanced and all considered only one response. In the present work, using two estimation algorithms was a complementary investigation and a proper comparison would require simulations from known parameters and models. Taking these limitations into account, some differences in parameter estimates were observed. Among the two software, the absorption rate estimates were those that differed the most, but they were also the most difficult to evaluate. The lack of sampling times early after dose administration led to poorly estimable absorption rate constants. Yet, we chose not to fix these parameters for the structural model building because this was the first modelling of S33138 and S35424 pharmacokinetics performed in patients and this population can greatly differ from healthy volunteers. Also, we were comforted in our decision by the consistency in the BIC and other parameter estimates. Difficulties were met to obtain SEs with FOCE-I in NONMEM despite using several sets of different initial conditions, however we did not systematically investigate initial conditions (this could be done using hypercube sampling in the multi-dimensional space of parameter). Non-parametric bootstrap can be performed to obtain an empirical distribution of the parameter uncertainty but in the present work we have performed permutation based Wald tests to correct for the departure from asymptotic conditions. Both procedures are computationally intensive and are facilitated by the use of closed form solutions which provides an important gain of time. While NONMEM easily handles steady state concentrations with model encoded in ODE, MLXTRAN in MONOLIX version 2.4 requires the user to add dummy lines of dose to mimic the path to steady state. The difference in computation time using this coding approach would thus increase with the number of subjects and the time to reach steady state.

CYP2D6 is an important catalyst of the oxidation of various antipsychotic agents: chlorpromazine, thioridazine, risperidone and haloperidol [START_REF] Murray | Role of CYP pharmacogenetics and drug-drug interactions in the efficacy and safety of atypical and other antipsychotic agents[END_REF]. Recently, significantly higher risperidone and 9-hydroxyrisperidone (its active metabolite) plasma concentrations have been observed in Korean CYP2D6 PM [START_REF] Kang | Effects of CYP2D6 and CYP3A5 genotypes on the plasma concentrations of risperidone and 9-hydroxyrisperidone in Korean schizophrenic patients[END_REF]. Here, both the metabolite and the parent drug showed a higher plasma exposure (AUC) in CYP2D6 PM patients compared to EM patients. This indicates that the back-transformation of the metabolite to the parent affects the disposition of the parent drug although concerning only 6% of the dose (30). Hamrén et al. have also described the impact of metabolite decreased elimination on parent drug concentrations in the presence of a reversible metabolism; in patients with impaired renal function, the decreased renal elimination of the tesaglitar metabolite (acyl glucuronide) causes an increase of the metabolite in plasma, which leads to increased amounts of metabolite undergoing back-transformation and subsequently an accumulation of tesaglitazar [START_REF] Hamrén | Mechanistic modelling of tesaglitazar pharmacokinetic data in subjects with various degrees of renal function-evidence of interconversion[END_REF]. In this respect, despite some limitations of the model (i.e. its apparent difficulty in capturing the lower concentration profiles in patients), the present work opens some perspectives for the ongoing development of the antipsychotic. The occurrence of CYP2D6 PM is higher in Caucasians (5-10%) than in East Asians (about 1%) [START_REF] Bertilsson | Pronounced differences between native Chinese and Swedish populations in the polymorphic hydroxylations of debrisoquin and S-mephenytoin[END_REF] and the CYP2D6 activity is lower in Chinese EM compared to Caucasians [START_REF] Bertilsson | Pronounced differences between native Chinese and Swedish populations in the polymorphic hydroxylations of debrisoquin and S-mephenytoin[END_REF] due to the CYP2D6*10 allele. Thus, additional investigations should be useful to assess the impact of CYP2D6 polymorphisms in populations other than Caucasians. As the model does not correctly predict the variability in healthy volunteers, it may not be used for comparison of Caucasian and Asian healthy volunteers.

However, it may be used as a model for Caucasian patients. Indeed, despite the lack of intravenous data in patients and identifiability problems, the present model give sensibles predictions of the concentrations of the parent drug and the metabolite in patients and of the effect of CYP2D6 on plasma exposure.

CONCLUSION

Using both algorithms and coding, we developed a complex joint pharmacokinetic model of interest for further developments of an antipsychotic and its metabolite. This model enabled us to evidence the CYP2D6 polymorphisms influence on the elimination of the active metabolite and adequately predict its impact on the parent drug levels through a back-transformation mechanism. V = V p = V m ω: between-subject standard deviation γ: within-subject standard deviation a: additive coefficient in nmol.L -1 b: multiplicative coefficient a p was fixed to 0 * Estimates that rely on the assumption that

LEGEND TO FIGURES

F p * K ap K am V p (=V m )* CL po * CL pm * CL mo * CL mp * (h -1 ) (h -1 ) (L) (L.h -1 ) (L.h -1 ) (L.h -1 ) (L.h -1 ) Without interconversion FOCE-I ODE - -99.0 NE - - 25
V m = V p C m = f DK a ((E 1 -K a ) + k pm F p /(1 -F p )) (V m /(1 -F p ))(λ 1 -K a )(λ 2 -K a ) e -Kat + f DK a ((E 1 -λ 1 ) + k pm F p /(1 -F p )) (V m /(1 -F p ))(K a -λ 1 )(λ 2 -λ 1 ) e -λ 1 t + f DK a ((E 1 -λ 2 ) + k pm F p /(1 -F p )) (V m /(1 -F p ))(K a -λ 2 )(λ 1 -λ 2 ) e -λ 2 t (6)
In both equations ( 5) and ( 6), E 1 is the parent drug total rate constant of elimination (= k po + k pm ), E 2 is the metabolite total rate constant of elimination (= k mo + k mp ), and λ 1 and λ 2 are the initial and terminal slopes of elimination, respectively defined in equations ( 7) and [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF]. 7)

λ 1 = (E 1 + E 2 ) + (E 1 + E 2 ) 2 -4(E 1 E 2 -k mp k pm ) 2 (
λ 2 = (E 1 + E 2 ) -(E 1 + E 2 ) 2 -4(E 1 E 2 -k mp k pm ) 2 (8)
Parameter identifiability

From the model slopes, the following parameters may be identified: K a , λ 1 and λ 2 .

From the model intercepts, we can identify the following equations:

V = (E 2 -λ 1 ) + k mp (1 -F p )/F p V p /f F p (9)
W = (E 2 -λ 2 ) + k mp (1 -F p )/F p V p /f F p (10)
Y = (E 1 -λ 1 ) + k pm F p /(1 -F p ) V m /f (1 -F p ) (11)
Z = (E 1 -λ 2 ) + k pm F p /(1 -F p ) V m /f (1 -F p ) (12)
From ( 7) and ( 8) we can write:

λ 1 + λ 2 = E 1 + E 2 (13)
λ 1 λ 2 = E 1 E 2 -k pm (F p /(1 -F p ))k mp ((1 -F p )/F p ) (14)
This yields a system of 6 equations, where V, W, Y and Z are four reals like λ 1 and λ 2 . We can estimate the following parameters:

V p /f F p , V m /f (1 -F p ), E 1 = k po + k pm , E 2 = k mo + k mp ,
k pm F p /(1 -F p ) and k mp (1 -F p )/F p , with the following solutions:

E 1 = (λ 1 + λ 2 -λ 1 λ 2 /T 2 -T 1 )/(1 -T 1 /T 2 ) E 2 = λ 1 + λ 2 -E 1 k mp (1 -F p )/F p = (λ 2 -λ 1 (W/V ))/(1 -W/V ) -E 2 k pm F p /(1 -F p ) = (λ 2 -λ 1 (Z/Y ))/(1 -Z/Y ) -E 1 V p /f F p = (E 2 -λ 1 ) + k mp (1 -F p )/F p V V m /f (1 -F p ) = (E 1 -λ 1 ) + k pm F p /(1 -F p ) Y (15)
where

T 1 = ( λ 2 W -λ 1 V )/( 1 W -1 V ) and T 2 = ( λ 2 Z -λ 1 Y )/( 1 Z -1 Y ).
Alternatively, the following parameterization may be used instead: [START_REF] Brendel | Are population pharmacokinetic and/or pharmacodynamic models adequately evaluated? A survey of the literature from 2002 to 2004[END_REF]

V p /f F p V m /f (1 -F p ) CL ptot /f F p = E 1 × V p /f F p CL mtot /f (1 -F p ) = E 2 × V m /f (1 -F p ) CL pm /V m = k pm V p /V m = k pm F p /(1 -F p ) × V p /f F p × f (1 -F p )/V m CL mp /V p = k mp V m /V p = k mp (1 -F p )/F p × V m /f (1 -F p ) × f F p /V p

  phase II study data set combining the observed profiles at both occasions contained 101 patients. Thirty five patients had a dose of 5 mg, 31 a dose of 10 mg and 35 a dose of 20 mg. The average [range] age in the population was 40 years [22-64], the average weight 69.0 kg [43.6-120.0]

For the two models

  without first-pass effect, the FOCE-I estimates were always included in the confidence intervals of the SAEM estimates. Whereas for the model with first-pass effect but no dose apportionment, the FOCE-I estimates for the volumes and all clearance but CL mp * did not fall within the confidence intervals of the SAEM estimates. For the model with dose apportionment and similar absorption rates, the confidence intervals of the FOCE-I estimates in ODE always included the FOCE-I estimates in CF and overlapped with the confidence intervals of the SAEM estimates. However with FOCE-I, the clearance of the parent by other pathways (CL po *) was greater than the clearance of transformation of the parent into the metabolite (CL pm *) while with SAEM CL po *<CL pm *. Also with FOCE-I, the clearance of the metabolite by other pathways (CL mo *) was similar to the clearance of transformation of the metabolite into the parent (CL mp *) while with SAEM CL mo *>CL mp *.Table III displays the between-subject standard deviation estimates of the four structural models as well as the model with similar absorption rates. The between-subject standard deviation estimates showed greater discrepancies using ODE or CF for the same model and estimation algorithm. Indeed, the median and maximal relative difference on clearance standard deviation estimates were respectively 0 and 150% with FOCE-I in NONMEM and 1.8 and 44% with SAEM in MONOLIX. As for the similarity based on the confidence intervals of the FOCE-I and/or SAEM estimates, the pattern was similar to that for the fixed effect. Of note with FOCE-I, ω CLpm * and ω CLmp * estimates were very close to zero, while zero was either included in, or very close to the boundaries of the corresponding confidence intervals of the corresponding SAEM estimates.The model selection proceeded similarly using either FOCE-I in NONMEM or SAEM in MONOLIX and ODE or CF, with the selected model including a back-transformation mechanism and a first-pass effect with a F p parameter for the dose apportionment. Both absorption/formation rate constants were difficult to estimate across models due to the combination of a 10 times ratio with the elimination rates and limited plasma data early after dose administration. Thus, we investigated a reduced model with similar absorption/formation rate constants for the parent drug and the metabolite. The BIC decreased by 30 units so that we kept this model for subsequent analyses. The RSE on the absorption rate constant in this model were about 15% with FOCE-I in NONMEM in ODE and 46 and 54% with SAEM in ODE and CF, respectively. Given the shorter computation time and the possibility to obtain standard error estimates, the MONOLIX software with the model encoded in closed form solution was used in the following.

Figure 6

 6 Figure 6 displays the npde from the model including the genetic covariate versus time for the parent drug and the metabolite on both occasions in the original data set. For the parent, the npde at 3, 6 and 24h were evenly distributed around 0 with the colored area formed by the 90% interval almost confounded with the dashed lines that represent the 5 th and 95 th percentiles of the normal

Figure 1 .

 1 Figure 1. Chemical structure of the parent drug: S33138 (a) and its active metabolite: S35424 (b).

Figure 2 .

 2 Figure 2. Schematic representation of the four tested structural models. f=fraction of dose after absorption, D=the dose, F p =fraction of parent reaching systemic circulation after absorption, K ap =absorption constant for the parent, K am =absorption constant for the metabolite, V p =volume of distribution for the parent, V m =volume of distribution for the metabolite, CL po =clearance of the parent by other pathways, CL pm =clearance of the parent into the metabolite, CL mo =clearance of the metabolite by other pathways, CL mp =clearance of back-transformation of the metabolite into the parent drug.

Figure 3 .

 3 Figure 3. Spaghetti plot of the observed concentrations of parent drug and metabolite versus time collected on two occasions, four and eight weeks after the treatment onset (W4 and W8) on a normal (a, b) and semi-log scale (c, d). The solid lines represent profiles from patients with a dose of 5 mg, dashed lines profiles from patients with a dose of 10 mg and dotted lines profiles from patients with a dose of 20 mg.

Figure 4 .

 4 Figure 4. Confidence interval visual predictive check plots of the selected structural and variability model including the dose effect, on a semi-log scale. The 90% confidence interval around the 5 th , 50 th and 95 th prediction percentiles from 250 simulated data sets are overlaid on the observed 5 th , 50 th and 95 th percentiles for the parent drug (top) and the metabolite (bottom) at W4 and W8 for a dose of 5 mg (a), 10 mg (b) and 20 mg (c).

Figure 5 .

 5 Figure 5. Boxplots of the parent drug (a) and metabolite (b) concentrations in the study, atthe three dose levels. The 10 th percentile corresponds to the lower whisker, the 25 th percentile to the lower hinge, the median to the thick bar, the 75 th percentile to the upper hinge and the 90 th percentile to the upper whisker.

Figure 6 .

 6 Figure 6. Normalized prediction distribution errors versus time from the final covariate model using the model building data set for the parent drug (top) and the metabolite (bottom) at W4 and W8. The corresponding 90% interval and the median are overlaid on the plot, and the dashed lines represent the 90% interval and median of the normal distribution.

Figure 7 .

 7 Figure 7. Normalized prediction distribution errors versus time from the final covariate model using the external evaluation data set for the parent drug (top) and the metabolite (bottom). The corresponding 90% interval and the median are overlaid on the plot, and the dashed lines represent the 90% interval and median of the normal distribution.

Figure 8 .

 8 Figure 8. Classic visual predictive check plots of the final covariate model using the external evaluation data set for the parent drug (top) and the metabolite (bottom), on a semi-log scale. The plain lines in grey represent the 5 th , 50 th and 95 th predicted percentiles and the dark dots jointed by dashed lines represent the 5 th , 50 th and 95 th observed percentiles for the parent drug (top) and the metabolite (bottom) for a dose of 10 mg (a), 20 mg (b) and 30 mg (c).

  genetic effect analysis indicated that CL mo * was decreased by 34% (p-value=0.015, Wald test by permutation) in CYP2D6 PM patients The population parameters of the model including the genetic covariate and their RSE are given in TableIV, for the 99 patients with available genotyping for the CYP2D6 polymorphisms. Based on these estimates for a CYP2D6 EM patients at a dose of 10 mg, the following estimates could be derived for the apparent parameters (with between-subject, and when applicable, within subject standard deviation in brackets) V p /F p =22.34 (0.42) L, V m /(1-F p )=145.27 (0.42) L, CL ptot /F p =3.19 (1.58 and 0.83) L.h -1 , CL mtot /(1-F p )=4.40 (0.65) L.h -1 , CL pm /V m =0.11 (0.28) h -1 and CL mp /V p =0.004 (0.28) h -1 . Figure5represents the distribution of the concentrations at each dose level between PM and EM. The rise in metabolite concentration levels in CYP2D6 poor metabolizers compared to extensive metabolizers appeared clearly in all dose groups while for the parent drug, the impact was only noticeable at the 20 mg dose. The between-subject standard deviation for CL mo * was only slightly decreased with the incorporation of the covariate. However, it is noteworthy that there were only 12 CYP2D6 PM patients in this study. Following Özdemir definition (26), R GC was computed only for CL po * because it was the only parameter with non-null within-subject variance. In this study, the R GC

for CL po * was equal to 0 as γ 2

CLpo * was superior to ω 2 CLpo * and indeed no genetic effect was found on this parameter. No effect of the CYP2C19 polymorphisms was found, probably due to the small number of PM.

Table I .

 I Bayesian Information Criteria (BIC), residual error model coefficient estimates and the computing time for the four structural models and the model with similar absorption rates, encoded in ordinary differential equations (ODE) and closed form (CF) using both the FOCE-I and the SAEM algorithms implemented in NONMEM and MONOLIX respectively, on data at W4

	BIC a p b p a m b m time	Without interconversion FOCE-I ODE 10 357 5.82 0.35 89.2 0.05 0.31	CF 10 368 4.00 0.36 96.4 0.05 0.14	SAEM ODE 10 352 27.60 0.30 105.0 0.05 0.98	CF 10 359 28.80 0.30 89.2 0.05 0.01	With interconversion FOCE-I ODE 10 241 0.34 0.34 76.2 0.05 0.50	Without first-pass effect CF 10 260 1.64 0.35 85.8 0.05 0.16	SAEM ODE 10 253 0.06 0.34 74.1 0.06 1.37	CF 10 245 0.02 0.34 70.0 0.05 0.01	With first-pass effect FOCE-I ODE 10 196 2.84 0.34 71.6 0.05 0.60	Without dose apportionment CF 10 206 2.97 0.35 74.9 0.05 0.27	SAEM ODE 10 212 0.08 0.34 74.8 0.05 1.40	CF 10 215 0.37 0.34 75.3 0.05 0.01	With dose apportionment FOCE-I ODE 10 186 2.47 0.34 55.90 0.05 53.8	Different absorption rates CF 10 196 2.44 0.35 53.0 0.05 0.37	SAEM ODE 10 176 0.57 0.32 58.3 0.05 2.42	CF 10 196 0.27 0.33 54.0 0.05 0.01	With dose apportionment FOCE-I ODE 10 149 3.65 0.31 51.30 0.05 35.4	Similar absorption rates CF 10 166 3.17 0.32 51.20 0.05 0.45	SAEM ODE 10 141 -0.01 0.32 57.90 0.05 2.37	CF 10 152 0.09 0.32 52.20 0.05 0.01	a: additive coefficient in nmol.L -1	b: multiplicative coefficient	time: computing time in hour

Table II .

 II Fixed effect estimates and their RSE (%, in grey) for the four structural models and the model with similar absorption rates, encoded in ordinary differential equations (ODE) and closed form (CF) using both the FOCE-I and the SAEM algorithms implemented in NONMEM and MONOLIX respectively, on data at W4

Table III .

 III Between subject standard deviation and their RSE (%, in grey) for the four structural models and the model with similar absorption rates, encoded in ordinary differential equations (ODE) and closed form (CF) using both the FOCE-I and the SAEM algorithms implemented in NONMEM and MONOLIX respectively, on data at W4

	Vm ) * ω CLpo * ω CLpm * ω CLmo * ω CLmp *	NE 2.04 NE 0.30 NE 0.51 NE --	NE 0.19 NE 2.66 NE 0.61 NE --NE 1.70 NE 0.29 NE 0.51 NE --	NE 0.28 NE 2.56 NE 0.59 NE --14 1.46 12 0.30 15 0.44 10 --	4 0.57 27 2.06 6 0.46 7 --13 1.60 18 0.37 12 0.52 8 --	4 0.31 39 2.40 6 0.54 7 --NE 0.95 NE 0.20 NE 0.60 NE 0.05 NE	NE 1.28 NE 2.70 NE 0.38 NE 0.17 NE NE 0.88 NE 0.21 NE 0.60 NE 0.00 NE	NE 1.37 NE 2.67 NE 0.38 NE 0.17 NE 17 1.09 12 0.21 27 0.62 9 0.31 40	4 1.00 22 2.95 6 0.40 11 0.20 9 12 0.98 13 0.24 22 0.59 9 0.28 50	5 0.99 21 2.91 7 0.42 10 0.18 9 NE 0.55 NE 0.02 NE 0.90 NE 0.01 NE	NE 2.55 NE 1.34 NE 0.18 NE 0.19 NE NE 0.55 NE 0.05 NE 0.87 NE 0.01 NE	NE 2.53 NE 1.31 NE 0.18 NE 0.19 NE 12 1.12 13 0.20 30 0.55 9 0.36 56	4 0.70 24 2.36 7 0.48 8 0.11 15 14 0.90 12 0.23 28 0.57 9 0.29 62	4 1.02 20 2.39 8 0.42 10 0.16 12 NE 0.54 NE 0.00 NE 0.91 NE 0.00 NE	NE 2.54 NE 1.21 14 0.18 NE 0.18 NE NE 0.54 NE 0.00 NE 0.86 NE 0.00 NE	NE 2.45 NE 1.24 NE 0.20 NE 0.18 NE 12 0.78 10 0.20 46 0.64 9 0.30 51	5 1.39 15 1.76 9 0.34 12 0.13 12 12 0.75 11 0.15 89 0.60 9 0.41 30	5 1.44 15 1.80 10 0.35 12 0.15 12 37 0.53 76 0.06 3579 0.91 54 0.01 117	10 2.58 34 1.16 41 0.18 66 0.18 21 NE 0.54 NE 0.05 NE 0.85 NE 0.01 NE	NE 2.48 NE 1.20 NE 0.19 NE 0.17 NE 12 0.10 10 0.20 36 0.54 10 0.34 43	4 1.15 18 2.06 8 0.39 10 0.14 11 12 0.12 12 0.14 61 0.53 9 0.49 28	4 0.80 23 2.33 7 0.46 8 0.13 14
	Vp (= ω	0.21	.2 0.20	25.0 0.28	23.0 0.31	23.7 0.27	20.9 0.25	20.9 0.23	20.9 0.31	20.8 0.28	17.4 0.28	17.3 0.30	20.8 0.27	20.8 0.31	17.0 0.33	17.0 0.36	18.4 0.30	18.3 0.34	16.5 0.35	16.6 0.31	18.9 0.29	19.4
	ω
	ω Fp * ω Kap ω Kam	Without interconversion FOCE-I ODE --0.05 NE --	CF --0.06 NE --	CF --99.0 NE --SAEM ODE --2.18 47 --	SAEM ODE --46.7 133 --CF --1.60 61 --	CF --29.2 120 --With interconversion FOCE-I ODE --3.35 NE --	With interconversion FOCE-I ODE --74.9 NE --Without first-pass effect CF --0.04 NE --	Without first-pass effect CF --99.0 NE --SAEM ODE --2.66 39 --	SAEM ODE --65.6 132 --CF --2.94 46 --	CF --101.0 171 --With first-pass effect FOCE-I ODE --0.01 NE 0.32 NE	With first-pass effect FOCE-I ODE --1.5 NE 0.3 NE Without dose apportionment CF --0.01 NE 0.29 NE	Without dose apportionment CF --1.5 NE 0.3 NE SAEM ODE --0.66 70 0.66 57	SAEM ODE --5.4 52 0.5 44 CF --0.68 96 0.76 78	CF --11.4 77 0.7 71 With dose apportionment FOCE-I ODE 0.00 NE 0.00 NE 1.37 NE	With dose apportionment FOCE-I ODE 0.82 NE 1.87 NE 1.88 NE Different absorption rates CF 0.00 NE 0.00 NE 1.42 NE	Different absorption rates CF 0.82 NE 1.92 NE 2.05 NE SAEM ODE 0.44 46 2.13 32 1.51 57	SAEM ODE 0.86 2 20.2 84 5.9 87 CF 0.12 21 2.00 39 4.89 21	CF 0.83 3 22.0 101 7.6 121 -With dose apportionment FOCE-I ODE 0.01 235 1.12 52 -	With dose apportionment FOCE-I ODE 0.81 2 1.8 40 ---Similar absorption rates CF 0.01 NE 1.15 NE -	Similar absorption rates CF 0.81 NE 2.0 NE ----34 1.66 22 SAEM ODE 0.07	SAEM ODE 0.87 2 7.5 46 ----CF 0.07 37 1.96 22	CF 0.85 2 9.3 54 --NE: not estimated	NE: not estimated ω: between-subject standard deviation	* Estimates that rely on the assumption that V

m = V p * Estimates that rely on the assumption that V m = V p

Table IV .

 IV Population pharmacokinetic parameters of the parent drug and its active metabolite for the model without and with the genetic covariate: estimates and relative standard errors (in grey, %) using the SAEM algorithm in MONOLIX and encoding in closed form

			Without	With
		the genetic covariate	the genetic covariate
			(N=101)	(N=99)
	f	1.00	-	1.00	-
	β f,D * (nmol -1 )	-0.02	28	-0.02	24
	F p *	0.85	2	0.87	2
	β Fp,D * (nmol -1 )	-0.04	29	-0.06	27
	K a (h -1 )	8.06	-	8.06	-
	V * (L)	19.2	4	19.4	5
	CL po * (L.h -1 )	0.84	11	0.67	11
	CL pm * (L.h -1 )	1.94	5	2.09	5
	CL mo * (L.h -1 )	0.46	6	0.50	7
	β CLmo,CY P 2D6 * (log(L.h -1 ))	-	-	-0.42	40
	CL mp * (L.h -1 )	0.09	12	0.09	12
	ω f *	0.25	13	0.27	12
	ω Ka *	1.50	21	1.41	25
	ω V *	0.19	18	0.26	13
	ω CLpo *	0.52	29	0.46	38
	ω CLmo *	0.51	8	0.50	8
	γ CLpo *	0.79	12	0.82	12
	b p	0.31	3	0.30	3
	a m (nmol.L -1 )	66.9	14	66.5	14
	b m	0.06	9	0.06	9

ACKNOWLEDGMENTS

The authors would like to thank Pr. N. Holford for his clever inputs on the choice of the parameter's notation to facilitate the understanding of the different structural models. We would also like to thank the MONOLIX development team for making the MLXTRAN tool available for this work.

During this work, Julie Bertrand was supported by a grant from the Institut de Recherches Internationales Servier (France).

APPENDIX

In all the equations below, C p is the parent drug concentration in plasma and C m is the active metabolite concentration in plasma following a single oral administration of a dose D of the parent drug.

The ordinary differential equation system and the corresponding closed form solution correspond to the last model on the right hand side of Figure 2, with similar absorption rates

Ordinary differential equation system

In the ordinary differential equation system (4), f is the fraction of dose after absorption, D is the dose, F p is the fraction of parent reaching systemic circulation after absorption, K a is the absorption constant for the parent and the metabolite, V p is the volume of distribution for the parent, V m is the volume of distribution for the metabolite, k po is the parent rate constant of elimination by other pathways (= CL po /V p ), k pm is the parent rate constant of transformation into the metabolite (= CL pm /V p ), k mo is the metabolite rate constant of elimination by other pathways (= CL mo /V m ), and k mp is the metabolite rate constant of back-transformation into the parent

Closed form solutions

(5)