{'fraction_non_alphanumeric': 0.03961176324194284, 'fraction_numerical': 0.019204446472271767, 'mean_word_length': 4.366125083463165, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compartmental analysis of dynamic nuclear medicine data: regularization procedure and application to physiology

Compartmental models based on tracer mass balance are extensively used in clinical and pre-clinical nuclear medicine in order to obtain quantitative information on tracer metabolism in the biological tissue. This paper is the second of a series of two that deal with the problem of tracer coefficient estimation via compartmental modelling in an inverse problem framework. While the previous work was devoted to the discussion of identifiability issues for 2, 3 and n-dimension compartmental systems, here we discuss the problem of numerically determining the tracer coefficients by means of a general regularized Multivariate Gauss Newton scheme. In this paper, applications concerning cerebral, hepatic and renal functions are considered, involving experimental measurements on FDG-PET data on different set of murine models.


Introduction
Nuclear medicine imaging is a class of functional imaging modality that utilizes radioactive tracers to investigate specific physiological processes. Such tracers are in general short-lived isotopes that are injected in the subject's blood and linked to chemical compounds whose metabolism is highly significant to understand the function or malfunction of an organ. Positron Emission Tomography (PET) [1] is the most modern nuclear medicine technique, utilizing isotopes produced in a cyclotron and providing dynamical images of its metabolism-based accumulation in the tissues.
Applications of PET in the clinical workflow depend on the kind of tracer employed and on the kind of metabolism that such tracer is able to involve: in this paper we will make use of [ 18 F]fluoro-2-deoxy-D-glucose-PET (FDG-PET) data. FDG is largely used for PET, mainly in the case of oncological applications [2,3,4,5].
This paper describes a very general numerical scheme for the reduction of different compartment models of FDG metabolism, based on a regularized multivariate Gauss Newton algorithm [6]. As many other approaches to compartmental analysis [7,8,9,10,11,12,13], also this method realizes numerical optimization but in a peculiarly effective way: the matrix differentiation step required at some stage of the analysis is here performed analytically, thus avoiding time consuming numerical differentiation. We test the reliability of the approach in the cases of a synthetic dataset and of three sets of experimental measurements concerned with the FDG metabolism in the brain, in the liver, and in the kidneys, provided by a micro-PET scanner for small animals.
The plan of the paper is as follows. Section 2 briefly recall the general ncompartment model for FDG metabolism. Section 3 introduces the numerical method for model reduction. Section 4 describes the FDG-PET models for the cerebral, hepatic and renal physiology. Section 5 provides numerical and experimental validation of the approach. Our conclusions are offered in Section 6.
In the document, R + , R − , R * + , R * − respectively denote the set of non-negative real numbers, the set of non-positive real numbers, the set of positive real numbers, the set of negative real numbers. For non-negative integers p and q, p, q denotes the set of integer larger than or equal to p and smaller than or equal to q. For a positive integer n, the canonical basis of R n is denoted by (e p ) p∈ 1,n and M n (R) denotes the algebra of n × n matrices with coefficients in R. For Banach or Fréchet spaces E and F , the vector space of bounded operators from E to F is denoted by L(E, F ).

The n-compartment system
We introduce the more general case of an n-compartment system and then describe the regularized Multivariate Gauss Newton algorithm. Simulation results and real data results will be provided on the more reliable cases of two-compartment and threecompartment models. The identifiability results of [14] insures that in the cases we consider, we can have good expectations of robusteness and accuracy even though the initial guess is not well-chosen, as applications can prove. However, in more general cases, either it is known that no uniqueness holds or when no information on uniqueness is known, the robustness and accuracy of the algorithm might greatly rely on the choice of a good initial guess. We consider a system composed of n compartments as in Figure 1, typically an organ, where each compartment represents a functional behaviour of a tissue, or an anatomical district. Remark that although the generic term "compartment" is used, it does not necessary mean that each compartment is contained in a physical compartment distinguishable from the others. In fact, this even constitutes one of the main issue in the inverse problem of getting information on the system since each compartment can not be individually observed. A radioactive tracer is injected to a patient and for a compartment p ∈ 1, n , C p denotes the non-negative concentration function of the tracer in the compartment. The compartment p receives the radioactive tracer from the outside world at a constant non-negative rate k pe and a non-negative concentration function C pe and it excretes the tracer at a constant non-negative rate k ep in the outside world. The constant non-negative rate at which the compartment p receives the tracer from a compartment q = p is denoted k pq . The concentration functions (C pe ) p∈ 1,n 1 Figure 1. n-compartment system are supposed to be continuous. The evolution of the tracer concentrations in each compartment is then governed by the following linear system of ordinary differential equations with constant coefficientṡ with the initial conditions where, for p ∈ 1, n , k pp = − q =p k qp + k ep . That iṡ and the matrix M is given by PET-scan images allow to know the total amount of radioactive tracer in the organ. Each compartment contributes to the intensity of the image linearly with respect to the amount of tracer in the compartment, that is its volume times the tracer concentration. Hence, in practice, PET-scan data gives access tõ where α ∈ R * + n is a known constant vector. Thus, the inverse problem we consider is to recover the exchange rates K ∈ R n 2 +n , where with measures ofC. Note however that some exchange rates may be a priori known.

A Newton algorithm for the inverse problem
In the following of the document, for a positive integer n and K ∈ R n 2 +n , we denote bŷ K ∈ R n 2 the first n 2 components of K andǨ ∈ R n the last n components of K. For a positive integer n, we denote by M the following linear operator H p+n(q−1) , p, q ∈ 1, n , p = q, so that for all H ∈ R n 2 + , M(H) is the matrix defined in (3) for the parameters H. Consider now the linear operator Vec : M n (R) → R n 2 stacking the columns of a matrix A into a column vector Vec(A), that is The operator Vec • M is linear from R n 2 into itself, we denote by S ∈ M n 2 (R) its matrix. It is a n 2 × n 2 sparse matrix with 2n 2 − n non-zero elements which are equal to ±1. More precisely, n 2 − n entries of S are equal to 1 and n 2 entries of S are equal to −1, such as We denote by W the following linear operator where for allČ ∈ C 0 (R + , R) n , H ∈ R n and p ∈ 1, n , [W(Č)(H)] p = H pČp , so that for all vector of input concentrations functionsČ = (C pe ) p∈ 1,n ∈ C 0 (R + , R + ) n and H ∈ R n + , W(Č)(H) is the vector defined in (2) for the input concentrations functionsČ and the parameters H. We denote by C the following function where for allČ ∈ C 0 (R + , R) n and K ∈ R n 2 +n , C = C(Č)(K) ∈ C 1 (R + , R) n is the unique solution tȯ

Algorithm description
We consider a n-compartment system as in Figure 1 with known input concentration functionsČ = (C pe ) p∈ 1,n ∈ C 0 (R + , R + ) n and α ∈ R * + n . We recall that the problem we are interested in consists in recovering the exchange rates K ∈ Ω ⊂ R n 2 +n , where K is given by (5), from the knowledge ofC(K), whereC = C α (Č). We recall that where C = C(Č). The concentrations vector C(K) is the solution of the ordinary differential equationṡ . The solution C(K) to the system of ordinary differential equations is given by C : R n 2 +n → C 1 (R + , R) n can be easily seen to be differentiable and even analytic. In order to use a Newton algorithm, we need to compute its Fréchet derivative. More precisely, considering t ∈ R + , we will compute the gradient ofC t for all t ∈ R + , with respect to K, where for all K ∈ R n 2 +n ,C t (K) =C(K)(t). For all K ∈ R n 2 +n , the Fréchet derivative dC dK (K) ofC at K, bounded operator from R n 2 +n to C 1 (R + , R) is then given by For all t ∈ R + , the gradient ofC t is given by where ∇K denotes the gradient with respect toK and ∇Ǩ denotes the gradient with respect toǨ. SinceC t is linear with respect toǨ, we simply have for all K ∈ R n 2 +n where denotes the Khatri-Rao product. Compute now the gradient ofC t with respect toK. WritingC t as where for all K ∈ R n 2 +n and N ∈ M n (R) Hence, to compute the gradient ofC t with respect toK, we need to compute the Fréchet derivative of F t with respect to the second variable N ∈ M n (R). The Fréchet derivative of the exponential function can be written in this way Hence, for all K ∈ Ω and N, H ∈ M n (R) so that making the change of variables (t − τ ) → τ , we get In other words Writing the previous formula in a more convenient way for computations, we have where ⊗ denotes the Kronecker product and the linear operator Vec : M n (R) → R n 2 , stacking the column of a matrix into a column vector, is defined in (6). Consequently, from (10), we have for all K ∈ R n 2 +n and H ∈ R n 2 Recalling that S, defined in (7), denotes the matrix of the linear operator Vec • M : Hence, using (8), (9), (11) In other words, for all K, H ∈ R n 2 +n and t ∈ R + , we have Consider a known functionC meas ∈ C 1 (R + , R), the measures ofC(K exact ) where K exact are the real unknown exchange rates of the compartmental system. Let K 0 ∈ Ω be an initial guess, then the Newton algorithm consists in solving the linear equation with Then, increment the value K 0 , giving K 1 = K 0 + H 0 and iterate the process.

Implementation details
The equation (12) may have no solution, moreover, in real applications, one has only the measured dataC meas for a finite number of sampling time points t 1 , . . . , t m ∈ R + and the data may be noisy. The discretized Newton algorithm consists in solving the linear system (12) by Tikhonov regularization. The non-regularized discretized system is given by and by Y 0 the vector the non-regularized discretized system can be written as Solving the system (13) by Tikhonov regularization consists in finding the solution H 0 to where r is a regularization parameter. As previously described, the value K 0 is incremented, giving K 1 = K 0 + H 0 and the process is iterated. Note that the regularization parameter can be different at each iteration step. As already mentioned, FDG is largely used for PET, mainly in the case of oncological applications. In fact, FDG-PET is based on the higher glycolytic activity in tumor cells compared to healthy tissue as the cause of image contrast. This glucose-analog behaviour implies that FDG is transported into malignant cells, which therefore exhibit increased radioactivity. From a biochemical viewpoint, FDG may follow a two-destiny path: in the first path FDG molecules are phosphorylated by means of a 6-phosphate group and remain trapped within the cells while, in the second path, FDG is not metabolized and therefore remains free in the tissue. This behaviour implies that, from a compartmental perspective, the analysis of FDG-PET data relies on a model made of an Input Function (in general, but not exclusively, an arterial one) and two functional compartments describing the free tracer and the trapped, metabolized FDG-6P, respectively.

The FDG-PET models
In the two-compartment model of FDG-PET data (see Figure 2) the state variables are the tracer concentration in compartments C f and C m for the free and metabolized compartments, respectively, while the kinetic process in the system is initialized by the Input Function (IF) C b , representing the tracer concentration in blood. The four constant transmission coefficients between the communicating compartments are denoted, for sake of simplicity, as k 1 , k 2 , k 3 , k 4 . The model equations for the forward and inverse problem are easily obtained from the general form in eq (1)-(5).
Next subsections describe three compartmental models for the FDG related to different physiological systems.

The model of the brain
For the FDG model of the brain, we will use a two-compartment catenary compartmental model. Indeed, even if from a functional (neuronal) viewpoint the brain has an extremely complex behaviour, brain-glucose metabolism inside the brain can be thought as standard [15]. We are aware that recently has been shown that the brain exchange coefficients of the compartmental systems, can vary very much according to the spatial position in the organ. For this reason, the two-compartment catenary compartmental system used for describing the brain physiology is modelled, applied and reduced pixelwise [16,17,18]. Such compartmental models are known as parametric compartmental models or indirect parametric imaging.
We will not discuss such parametric models, and otherwise discuss the validity and robustness of our technique if applied to a two-compartment system, using brain data to validate it on real imaging FDG-PET data. For liver studies, the FDG-PET compartmental analysis is complex, due to the hepatocyte capability to de-phosphorilate and release glucose and FDG into the blood. Moreover, a tracer is supplied to the liver by both the hepatic artery and the portal vein (dual input). Therefore the compartmental model that we need to adopt, is a two-compartment catenary one, with two IFs. However, the portal vein is not visible in PET images; in order to overcome this limitation we coupled the capability of PET to provide FDG concentration curves in virtually all the organs with the anatomical ken that the portal vein enters the liver just after leaving the gut, as in [19]. Indeed, we use this information to reduce a compartmental scheme in which gut FDG metabolism is modelled as a standard two-compartment compartmental system, its tracer concentration is used to compute the portal vein IF, and therefore a model as in Figure 3 is reduced, to describe tracer kinetics in the hepatic system. Clinical activity asks to represent glucose consumption by malignant cells as Standardized Uptake Value (SUV) [20] whose estimation is hampered by the interference of several variables modulating tracer availability independently of tumor metabolism [21]. Among these factors, kidney function plays a relevant role as, differently from the tracked glucose, FDG is scarcely reabsorbed in renal tubule and is largely excreted in the urines [22]. Furthermore, variations in FDG concentration due to water reabsorption in renal tubules and increase of bladders volume during the FDG excretion process have to be taken in account, in order to describe the renal FDG metabolism: therefore we describe the renal system as a three-compartmental system. According to [8], we have some physiological constraints that simplify the compartmental scheme: k 7 is proportional to k 4 with a known proportional factor λ = 10 2 , i.e. k 4 = 10 2 k 7 ; k 5 is set to zero since the tracer is not reabsorbed in tubules and goes to bladder (which is the output k 7 of the compartmental scheme).

Numerical Results
An 'Albira' micro-PET system produced by Carestream Health is currently operational at the IRCCS San Martino IST, Genova. Experiments with mice are currently performed at this site, using different tracers, mainly for application to oncology and to FDG physiology in the renal and liver system [23,7,8,19]. In this section we describe the performance of our approach to compartmental analysis in the case of both synthetic and real 'Albira' data.

Synthetic data
We realized three different simulation experiments, one for each model described in this paper: a two-compartment scheme for the brain studies, a two-compartment with dual 1.00 0.20 0.05 0.80 MGN 1.05 ± 0.04 0.20 ± 0.02 0.05 ± 0.01 0.82 ± 0.03 LS 1.02 ± 0.08 0.18 ± 0.05 0.05 ± 0.05 0.80 ± 0.03 Table 1. Validation with a synthetic dataset: ground-truth (g.t.) values of the tracer coefficients compared with the reconstructions provided by multilinear fitting with the Levenberg-Marquardt Least-Squares algorithm (LS) and by our regularized Multivariate Gauss Newton (MGN) approach. V is set equal to 0.02.
input model for the liver (coupled with a standard two-compartment scheme for the gut) and a three compartment model for renal studies.
In order to produce the synthetic data we initially chose realistic values for V (the blood fraction that supplies the organ of interest) and for λ of the renal system; then we utilized ground truth (g.t.) values for all the tracer kinetic parameters (4 in the brain case, 4+4 in the gut-liver case and 5 in the kidneys case). With these selected values we solved every forward problems (1)-(3) in term of C. These solutions of the Cauchy problems are sampled in time, on a total interval [t 1 , t m ] that corresponds to the total acquisition time with 'Albira' and in correspondance of time points typical of experiments with the scanner in this application context. The Input Function C b has always been obtained by fitting with a gamma variate function a set of real measurements acquired from a healthy mouse in a very controlled experiment [24]. We computed the dataC as in (4), affected the data by Poisson noise, and applied our algorithm in order to reconstruct the exchange coefficients. The Tikhonov regularization parameter r was optimized at each iteration using the Generalized Cross Validation technique [25]. Comparison with the g.t. values for this parameters provides limits about the reliability of the model and of the inversion procedure.
The results of this test are given in Tables 1, 2 and 3, where mean and standard deviations are computed over 50 runs of the same problem with different (random) initialization vectors K (0) , in a Monte Carlo approach. Iterations are stopped with a discrepancy principle on data in combination with a tolerance on the step-size H of the Newton Method. The computational burden is 10 seconds for each Monte Carlo run (on a Intel i5 2.3GHz x 4). In Tables, we show also comparison with a standard Levemberg-Marquardt method for the least-squares minimization [26]. Comparison with respect to the ground-truth values and the values provided by the Levenberg-Marquardt method, and the small values for the reconstruction uncertainties clearly show the reliability of our approach.

Real FDG-PET data
For our experiment we consider a control group (n=10) and a group (of the same size) in which FDG injection was performed after one month of high dose metformin treatment  Table 3. Validation with a synthetic dataset: ground-truth (g.t.) values of the tracer coefficients compared with the reconstructions provided by multilinear fitting with the Levenberg-Marquardt Least-Squares algorithm (LS) and by our regularized Multivariate Gauss Newton (MGN) approach. V is set equal to 0.3, and the proportional factor λ is 10 2 (such that k 7 = 1 λ k 4 ).
(750 mg/Kg body weight daily). This drug reduces blood glucose concentration without causing hypoglycemia [27], mostly by decreasing intestinal glucose absorption and glucose delivery by the liver [23]. It follows that nuclear medicine experiments with such kind of animal models are perfect candidates to validate a compartmental model that has been designed in order to follow the FDG kinetics in the hepatic and renal system in a refined fashion [8,19].
To ensure a steady state of substrate and hormones governing glucose metabolism, the small animal was studied after six hours fasting. Mouse was weighted and anaesthesia was induced by intra-peritoneal administration of ketamine/xylazine (100 and 10 mg/kg, respectively). Serum glucose level was tested and animals were positioned on the bed of the scanner, whose two-ring configuration permits to cover the whole animal body in a single bed position. The mouse was injected with a dose of 3-4 MBq of FDG through a tail vein soon after the start of a list mode acquisition lasting 50 minutes. Acquisition was performed using the following framing rate: 10 × 15 sec, 5 × 30 sec, 2 × 150 sec, 6 × 300 sec, 1 × 600 sec. PET data were reconstructed using a Maximum Likelihood Expectation Maximization method (MLEM) [28]. Thereafter, each image dataset was reviewed by an experienced observer who recognized five Regions Of Interest (ROIs) encompassing left ventricle, brain, gut and liver and kidneys respectively (as in Figure 5). The ROIs over the left ventricle allowed us to compute the IF (we are aware that the determination of IF is a challenging task in the case of mice. To accomplish it, for each animal model we have first viewed the tracer first pass in cine mode; then, in a frame where the left ventricle was particularly visible, we have drawn a ROI in 0.982 ± 0.063 0.375 ± 0.051 0.010 ± 0.005 1.082 ± 0.093 metformin 1.051 ± 0.916 0.264 ± 0.031 0.004 ± 0.003 0.820 ± 0.032 Table 4. Mean and standard deviations for liver kinetic parameters over the two sets of mice (control, metformin).
1.027 ± 0.109 0.193 ± 0.027 0.187 ± 0.021 0.415 ± 0.76 0.512 ± 0.043 metformin 1.154 ± 0.110 173 ± 0.030 0.003 ± 0.003 0.490 ± 0.059 0.497 ± 0.050 Table 6. Mean and standard deviations for liver kinetic parameters over the two sets of mice (control, metformin). k 7 = 0.042 ± 0.017 the aortic arc and maintained it for all time points). The other ROIs allowed us to estimate the input data for the compartmental method, i.e. estimates ofC(t) all experimental time points. As far as V is concerned, we utilized physiological sound values for each ROI, accordingly to [29,30]. The output of our algorithm for model reduction is described in Table 4, 5 and 6, which contains the mean values and standard deviations of the reconstructed tracer coefficients after 50 runs of the code for 50 random initialization vectors. In the hepatic case, only results for liver are shown, the one of actual physiological interest.
These results are physiologically plausible: for example, in all the three experiments, values of k 3 decrease in metformin-treated mice with respect to the control mice (the brain case is the less notable one), coherently with the property of metformin to decrease de-phosphorilation of FDG (and glucose). The relative behaviour of the coefficient values in the kidneys and liver case are coherent with the results of analogous experiments [8,19].

Conclusions
This paper, as a sequel of the previous work concerning the compartmental inverse problem, discusses the mathematical model for the compartmental analysis and provides a general technique for its resolution. In particular, we described the performance of a regularized Newton approach to the numerical solution of the inverse problem associated to a compartmental model for FDG metabolism in micro-PET experiments. The numerical method is able to determine with a notable level of robustness and reliability all the tracer coefficients embedded in the model and, in particular, it does not a priori forbid de-phosphorylation. From a mathematical perspective, a further, more fundamental advantage of our scheme with respect to the Levenberg-Marquardt approach, is that Newton methods search for zeroes of non-linear functionals, and therefore do not need to a priori select a topology in the data space, as in the case of least-squares approaches. Results on real data verify the reliability of the method in estimating physiologically sound values for the tracer coefficients, that are in agreement with literature.
A Matlab prototype implementing this approach is at disposal at together with a Graphical User Interface for user-friendly input/output processing.