Geodesic Uncertainty in Diffusion MRI

We study theoretical and operational issues of geodesic tractography, a geometric methodology for retrieving biologically plausible neural fibers in the brain from diffusion weighted magnetic resonance imaging. The premise is that true positives are geodesics in a suitably constructed metric space, but unlike traditional first order methods these are not a priori constrained to connect nongeneric points on subdimensional manifolds, such as the characteristics in traditional streamline methods. By virtue of the Hopf-Rinow theorem geodesic tractography furnishes a huge amount of redundancy, ensuring the a priori existence of at least one tentative fiber between any two points and permitting additional tractometric and data-extrinsic constraints for (fuzzy or crisp) classification of true and false positives. In our feasibility study we consider a hybrid paradigm that unifies existing ideas on tractography, combining deterministic and probabilistic elements in a way naturally supported by metric geometry. Particular attention is paid to an analytical prediction of geodesic deviation on numerically computed geodesics, a ‘tidal’ effect induced by small perturbations resulting from data noise. Taking these effects into account clarifies the inherent uncertainty of geodesics, while simultaneosuly offering a dimensionality reduction of the tractography problem.


INTRODUCTION
A hydrogen nucleus (proton) behaves like a tiny magnet when interacting with an external magnetic field. In this situation the Zeeman splitting effect discloses two quantum eigenstates of its nuclear magnetic moment, distinguishable by a relative energy gap. These states are known in the trade as "spin up", or "parallel", and "spin down", or "anti-parallel". The attributes refer to the relative alignment with the external field, but should be taken with a grain of salt for quantum systems of individual nuclei. In a typical macroscopic system, however, the Zeeman effect induces a relatively small yet measurable "classical" magnetization [Bloch (1946); Cowan (1997)], for which these attributes can be taken literally.
Magnetic resonance imaging (MRI) allows one to manipulate the magnetization of a hydrogen spin system in a way that enables non-invasive in-vivo acquisition of informative localized signals, i.e. "images". Although tuned in on water-bound hydrogen, because of its abundance in the human body, the recorded signals are affected by the physicochemical environment and therefore carry (indirect) anatomical and/or functional information depending on the measurement protocol. MRI is a versatile modality with endless options for image formation, cf. (Brown et al. (2014)).
The empirical feature of interest for our purpose is the anisotropic diffusivity of mobile (waterbound) hydrogen spins in the human brain 1 . The MRI protocol that allows one to infer this from measurable nuclear magnetic resonance and relaxation effects is known as diffusion weighted imaging (DWI) (Torrey, 1956;Moritani et al., 2009). Insight in local water diffusivity patterns is particularly interesting for studying fibrous tissues, since it is stipulated that anisotropic media impart non-random barriers for diffusion, inducing relatively long mean free paths along "preferred orientations" (candidate fibers). In particular, DWI is the only non-invasive in-vivo technique for probing white matter pathways (bundles of myelinated axons or nerve fibers, aka "tracts") in the brain. As such it plays a pivotal role in connectomics, the 21st century's grand challenge that aims to disclose structure and function of the human brain.
DWI produces massive, non-visual data, calling for the right questions to be posed in order to render such data insightful in congruity with our visual inclination. An important question pertains to tractography, the inverse problem that aims to infer white matter pathways from the observed local diffusivity patterns. Such tracts furnish the wiring between various grey matter areas containing nerve cell bodies, the brain's computational units so to speak. For this reason tractography also plays an essential role in clinical applications, such as neurosurgery, in which spatial maps of anatomical tracts and functional networks are relevant for effective patient management.
However, tractography has inherent limitations that should be taken into account as well, e.g. inter-and intra-rater variablity makes comparisons of tractograms a challenging endavour (Maier-Hein et al., 2017;Schilling et al., 2019). In a recent paper by Schilling et al. (Schilling et al., 2020) the authors argue that tractography is anatomically plausible if furnished by appropriate constraints, notably end point conditions and nogoes. This observation argues for a tractography paradigm that is compatible with end point conditions imposed by the user and that does not exclude the possibility of multiple connections, two properties not shared by classical tracing (streamline) methods based on diffusion tensor imaging (DTI). Moreover, no-go conditions should be at the discretion of an expert, and not a result from a priori methodological limitations. An example is the fractional anisotropy (FA) stopping criterion of classical streamline methods, which results in incomplete trajectories due to destructive interference of anisotropic diffusivities at fiber crossings or severely reduced effective anisotropies in grey matter.
The conceptual issues alluded to above are gracefully handled by geodesic tractography, the mathematical details of which are explained in section 2. It is important to realize that this is not a one-on-one replacement of streamline tractography, but a genuine generalization that (intentionally) raises a new problem: Any two points in the brain are connected by at least one geodesic ("geodesic completeness"). Thus geodesics should not be confused with (streamlines intended to reflect) biological fibers, but should instead be seen as elements of a highly redundant data representation. The premise is, of course, that these can be pruned using data intrinsic and external knowledge to (hopefully) produce the desired true positives without incurring false negatives. That is, the assumption is that biological tracts are geodesics, but not vice versa. Interestingly, one can show that classical streamlines emerge as geodesics when the two minor eigenvalues of the DTI tensor are set to zero. This Gedanken experiment connects the two tractography paradigms conceptually, but also shows (by an argument of continuity) that there exist "streamline-like" geodesics in brain regions characterized by high FA, with "streamline-unlike" infillings in isotropic regions.
For the reasons above we focus on geodesic tractography, and present some feasibility studies in Section 3. More specifically, we aim to establish a conceptual and operational framework geared towards clinical needs in a neurosurgical workflow for brain tumor patients, in collaboration with the Department of Neurosurgery of Elisabeth Tweesteden Hospital (ETZ), a highly specialised referral center for brain tumor surgery in Netherlands (Meesters et al., 2019;Meesters et al., 2021;Rutten et al., 2014). Our approach is motivated by conceptual as well as practical demands in this context, viz. 1) to overcome the rigidity and shortcomings of classical tracing methods in return for a more versatile one, and 2) to connect to clinically viable MRI protocols, notably DTI, while keeping tabs on the development of more sophisticated DWI models for future refinement.
We sketch the geodesic tractography rationale in all generality, which in principle requires so-called Finsler geometry, but work out details for its application to DTI only. This case is mathematically considerably less mind-boggling since it allows us to stay within the realm of well-established Riemannian geometry. Another practical issue is the effect of noise and errors in (automatic or manual) seeding, which can be naturally accounted for by considering the concept of geodesic deviation, furnishing geodesic tracks with volumetric counterparts in the form of tubular uncertainty neighbourhoods ("geodesic tubes") .

THEORY
Taking the empirical adequacy of DWI for tractography for granted, something must be optimal along biologically meaningful curves. This suggests a variational approach, in which tracts arise as extrema of some functional (Arnold, 1989;Lovelock and Rund, 1988;Rund, 1973;Sagan, 1992;Weinstock, 1974).
In geodesic tractography, the functional of interest is constructed on the basis of a norm or metric. The relevant geometry is known as Riemann-Finsler geometry, or Finsler geometry for short (Bao et al., 2000;Cartan, 1934;Shen and Shen, 2016;Szilasi et al., 2014). Inner product induced norms are quite special, and often used out of habit rather than necessity. Such norms fall within the much narrower realm of Riemannian geometry, an area much better understood and (therefore) exploited than the general framework (do Carmo, 1976;do Carmo, 1993;Cartan, 1963;Jost, 2011;Koenderink, 1990;Misner et al., 1973;Spivak, 1975). The essential constraint is that inner product norms are square roots of quadratic forms, characterized by a small number of degrees of freedom, viz. six per point in three spatial dimensions for specifying the independent parameters of the defining positive-definite symmetric Gram matrix. Such norms are especially suited for applications in which the observable of interest is compatible with the quadratic restriction, notably diffusion tensor imaging (DTI). General DWI models, however, do not have this limitation. Extension of the variational approach based on metric geometry beyond DTI would therefore render a genuine Finslerian approach natural and compulsory.
Although Riemannian geodesic tractography does not cover traditional "streamline tractography", the latter may be viewed as a singular limit that can be inferred from a sequence of Riemannian metrics. In turn, Finsler geodesic tractography may be related to an orientation-parametrized family of Riemannian metrics (Astola, 2010;Florack and Fuster, 2014;Florack et al., 2015;Dela Haije, 2017;Dela Haije et al., 2019). In this sense Riemannian geodesic tractography plays an intermediate role, being a generalisation of the original streamline case and a special instance of the Finslerian case. Its limitations are offset by the inspiration that can be drawn from the vast body of literature on Riemannian geometry, as well as by the particular clinical appeal of DTI.
Positive definite matrix fields play a central role in DWI, and in DTI in particular, and so do, a fortiori, positivity preserving operators. In a differential geometric approach, differential operators are of prime interest. Unfortunately, positivity is not manifest in "standard" calculus. A more suitable framework for differentiation and integration is provided by multiplicative calculus, which, in the non-commutative case of positive symmetric matrices, is a non-trivial extension in which positivity is manifest (Bashirov et al., 2008;Florack and van Assen, 2012;Gill and Johansen, 1990).
To make this article self-contained the most relevant mathematical concepts are loosely explained, de-emphasizing proofs, in subsections 2.1 (calculus of variations), 2.2 (multiplicative calculus) and 2.3 (Riemann-Finsler geometry). Literature pointers should help the reader interested in more detail.

Calculus of Variations
In the Lagrangian formalism of the calculus of variations the point of departure is an action functional, The integrand is called the Lagrangian. In many cases Lagrangians only depend on the m-dimensional independent variable x, the n-codimensional dependent variable ϕ, and its first order x-derivatives, ∇ϕ: The variational principle relates the physical manifestation of a Lagrangian system to critical points ϕ crit of the action functional, formally defined in terms of a vanishing functional derivative: Although the left hand side can be rigorously defined, it suffices to appreciate its intuitive meaning as the "rate of change" of the action functional S(ϕ) with respect to any small perturbation δϕ of the dependent variable ϕ relative to a fiducial instance ϕ crit . Eq. 3 yields a system of n differential equations known as the Euler-Lagrange system.
Our case of interest will be m 1, n 3. In this situation the domain of definition is a one-parameter interval, and one usually writes t ("time") or s ("length") instead of x, and L instead of L. Since in this case we have an ordinary rather than partial derivative operator, the ∇-symbol now represents d/dt or d/ds. Our codomain will be a 3-dimensional tangent space, representing the linear space of instantaneous "velocities'" within (any tangent space of) 3-dimensional space. For this reason one typically uses the spatial coordinate symbol x instead of ϕ. Its first order derivatives are then conveniently written either as _ x dx/dt (Newton's dot-notation) or as x′ dx/ ds (Lagrange's prime notation). In tractography the independent parameter serves as an arbitrary label for disambiguating points along a tract; the connotation with "time" or "length" is metaphorical.
Thus the action functional that concerns us here takes the form in which x x(t) represents the coordinate triple along a parametrized tract over an interval I ⊂ R. Stationarity entails which yields an Euler-Lagrange system in the form of three ordinary differential equations, one for each component _ x i , i 1, 2, 3, of the tangent vector _ x^_ x 1 e 1 + _ x 2 e 2 + _ x 3 e 3 relative to a fiducial local frame of basis vectors {e 1 , e 2 , e 3 }: The precise form of the Lagrangian is constructed on the basis of geometric considerations and will be specified later. It is important to realize that the spatial variability of the local frame (more precisely, its parallel transport 2 ) affects the form of the Euler-Lagrange equations.
For a thorough introduction to the calculus of variations, cf. Lovelock and Rund (1989). Rund elaborates on the profound implications of one-homogeneous Lagrangians Rund (1973). The book by Neuenschwander (2011) is particularly interesting in case the Lagrangian exhibits certain symmetries, explaining Noether's seminal work on this issue.

Multiplicative Calculus
Multiplicative calculus (Bashirov et al. (2008)) provides a natural framework in problems in which positive images or positive symmetric matrix fields and positivity preserving operators are of interest. In the commutative case it is a convenient, albeit arguably redundant framework, which can be restated in terms of standard (additive) calculus via a log-exp commuting diagram (Arsigny et al., 2006;Fillard et al., 2007;Florack and Astola, 2008;Pennec et al., 2006). Figure 1 illustrates this for the case involving a single independent variable x ∈ R, with the multiplicative derivative of a positive function f defined as The multiplicative antiderivative is then naturally defined as The left hand side notation is due to Gill (Gill and Johansen, 1990) in analogy with Leibiz "continuous-sum" symbol.
For a product antiderivative or indefinite product integral we have Definite product integrals are introduced via a spatial partitioning and limiting procedure, akin to the standard Riemann sum approximation in the additive case, b a f(x) dx lim The relationship between Eq. 9, and Eq. 11 is formalized by the fundamental theorem of multiplicative calculus: recall the well-known classical counterpart relating Eq. 10 and Numerical approximations of Eq. 11 are obtained by partitioning the integration domain and suppressing the limiting procedure.
It requires minor efforts to generalize foregoing results to the multivariate case. The non-commutative case, on the other hand, is highly nontrivial and not without ambiguity. The complication is apparent from the ambiguous ordering of non-commuting factors on the right hand side of Eq. 11. There are multiple options to establish an unambiguous definition, cf. Gantmacher (2001) and Slavík (2007)) We will briefly discuss some and subsequently pick the one suited for our purpose.
The log-exp framework may seem the most straightforward option, defining the multiplicative derivative of a positive symmetric matrix function according to Eq. 7, in which log and exp are replaced by their well-defined matrix counterparts. The chain rule for matrix functions, however, complicates matters. Declining from the assumption of commutativity we have, e.g.
for an arbitrary (square) matrix field g. For a positive symmetric matrix field f we furthermore have in which I is the identity matrix. Both identities can be readily proven with the help of the defining Taylor series expansions of the exp and ln functions, cf. Higham (2008). In the commutative case (only) Eqs 15 and 16 reproduce the familiar chain rule formulas FIGURE 1 | Commuting diagrams for (A) multiplicative differentiation, cf. (7) and (B) multiplicative integration, cf. (8) in terms of standard calculus counterparts via the log-exp route in the commutative case. The (upward left) operators p and are defined by virtue of the equivalent counterclockwise "classical" routes in these diagrams.
Frontiers in Computer Science | www.frontiersin.org October 2021 | Volume 3 | Article 718131 and well-known from classical calculus. In view of this one may avoid the cumbersome expression of Eq. 16 resulting from the log-exp trick of Eq. 7 by defining the multiplicative derivative in a different way, viz. by replacing the expression (ln f(x))′ on the right hand side in Eq. 7 formally by one of the alternative expressions in Eq. 18. As these two options (or any alternative construct obtained by sandwiching factorizations of f( (x)) are clearly different, the application context should determine one's bias. We will adhere to the following definition: It can be shown that, with this choice, Eq. 11 is applicable, provided we use a right-to-left ordering of factors on the right The reason for choosing Eqs 19 and 20 is that these provide the correct derivative/antiderivative for recasting a linear matrix differential equation f′ A f, with positivity constraint on f-which will be of interest later, notably in Eq. 48 on page 13-, into its natural multiplicative counterpart with built-in positivity preservation, viz. f* exp A. The solution then takes the trivial form of an antiderivative, f(x) exp(A(x)dx), according to 3 Eq. 20.
Another application of multiplicative calculus relevant for our purpose is the so-called log-Euclidean paradigm for symmetric positive-definite matrix functions 4 (Arsigny et al., 2006;Fillard et al., 2007;Pennec et al., 2006). We consider it in the context of multiscale representations (Florack and Astola, 2008). Figure 2 shows how to construct a consistent multiscale representation f t of a raw symmetric postive-definite matrix function f^f 0 , with t > 0 denoting (quadratic) scale, according to the following multiplicative convolution recipe: in which ϕ t denotes the normalized isotropic Gaussian of (quadratic) scale t and p acts component-wise: This scheme ensures closure and commutativity of blurring and inversion: If f t^e xp(ϕ t p ln f) and g t^e xp(ϕ t p ln g), with g f −1 , then g t f −1 t at all scales t > 0.

Riemann-Finsler Geometry
Riemann-Finsler geometry is a generalization of Riemannian geometry. Both are metric geometries, concerned with spaces endowed with a distance concept. The generalization entails the replacement of an inner product induced norm (square root of a quadratic form reflecting the Pythagorean rule) by a general norm. For this reason Riemann-Finsler geometry is sometimes referred to as 'Riemannian geometry without the quadratic restriction'. First studied by Finsler (1918), it has matured towards its present form mainly due to Cartan (1934). The mathematics is rather mind-boggling. For modern in-depth accounts cf. Bao et al. (2000), Shen and Shen (2016), and Antonelli et al. (1993); (Antonelli and Zastawniak, 1999) in relation to some applications in physics and biology. Riemann-Finsler geometry is also attracting considerable attention in the context of (quantum) gravity, see for example (Fuster and Pabst, 2016;Gibbons et al., 2007;Girelli et al., 2007;Pfeifer and Wohlfarth, 2012). For our purpose a norm is needed to introduce an abstract notion of "length" of a curve, one for which units of length are proportional to local mean free path lengths of diffusing water molecules in the underlying fibrous tissue, as observed via DWI. This implies that if local mean free paths are indeed longest along biological tracts, then such tracts would define "locally shortest" paths, or geodesics, along a Riemann-Finsler manifold. The attribute "local" means that any sufficiently small perturbation of a fiducial geodesic with fixed endpoints would incur an increase of length. This leaves room for multiple geodesics connecting the same pair of endpoints. The Hopf-Rinow theorem (Jost, 2011) ensures geodesic completeness, i.e. the existence of a "minimal geodesic", the length of which corresponds to (or, if you want, defines) the distance between its endpoints. In the case of multiple local minima between two endpoints, the global one is not necessarily the biologically most plausible one.
The Riemann-Finsler manifold is thus constructed in order to turn the tractography problem into a geodesic problem, which falls within the scope of the calculus of variations, recall section 2.1. The inherently positive nature of distance naturally calls for multiplicative calculus, recall section 2.2.
The most general functional expressing the Riemann-Finsler length of a tentative tract, arbitrarily parametrized as x x(t), has the following form, recall Eq. 4: in which the general Lagrangian L(t, x, _ x) has been replaced by the so-called Finsler norm F(x, _ x), expressing the Riemann-Finsler length of the tangent vector _ x at point x in a parameter invariant way. The latter not only means that it does not explicitly depend on the curve parameter t ∈ I, but also that it must be absolutely homogeneous of degree one, i.e. 5 for any tangent vector rescaling by λ∈ R.
One can show that, under mild conditions besides Eq. 24, there exists a positive-definite symmetric Riemann-Finsler metric tensor g ij (x, _ x), such that 6 The relation between F(x, _ x) and g ij (x, _ x) is, remarkably, one-to-one, with The Riemann-Finsler metric tensor is homogeneous of degree zero, Thus a Riemann-Finsler metric has an orientation dependency, but magnitude and "polarity" of the tangent vector _ x are immaterial. The quadratic restriction of Riemannian geometry entails a severe restriction on the Finsler norm, viz. g ij (x, _ x) g ij (x) independent of orientation. Equation 26 may be interpreted as a family of Riemannian metrics, one for each orientation. All members of this family clearly coincide in the Riemannian case, so that Eq. 25 reduces to its Riemannian counterpart, the Pythagorean form The Euler-Lagrange equations for Eq. 25 are known as the geodesic equations, and take the form The right hand side can be seen as a pseudoforce acting along the geodesic and as such merely reflects the arbitrariness of "time"-parametrization, contributing to an artificial acceleration along the trajectory. It has no effect on the relevant (parameter-invariant) geometry of the curve and can be effectively removed via (nonlinear) reparametrization. Any new parameter thus obtained is then coined an affine parameter, yielding "constant speed" geodesics characterized by F(x(t), _ x(t)) constant. The geodesic spray coefficients G i (x(t), _ x(t)) can be viewed as forces inducing an acceleration and bending of the trajectory that is completely determined by the inhomogeneous nature of the metric tensor, cf. Bao (Bao et al., 2000) or Shen and Shen Shen and Shen (2016) for technical details. In the Riemannian limit they assume a quadratic form in the velocities _ x. Figure 3 illustrates the rationale of constructing the action functional Eq. 23 as a data-induced curve length in the context of the variational principle for the Riemannian case applicable to DTI.
Apart from the geodesic equations, Eq. 29, we will be interested in "tidal effects', i.e. the relative separation acceleration 7 of neighbouring geodesics, aka geodesic deviation. This is relevant for understanding how stable a geodesic is with respect to small perturbations. By quantifying this one can analytically predict (to linear approximation) the behaviour of sufficiently narrow bundles of geodesics around any fiducial geodesic. Understanding geodesic deviation in a general metric space requires some fairly sophisticated mathematical machinery. It suffices here to convey the gist of it.
The geodesic deviation equations form a linear system of ordinary differential equations for a deviation vector ξ ξ(t) defined along a given geodesic curve x x(t). This vector can be seen as connecting points on two 'infinitesimally' neighbouring geodesics at the same parameter value t (one must, to begin with, first agree on an unambiguous affine parametrization). We have The δ-derivative is a so-called covariant modification of the dderivative from standard calculus, involving geometric 'correction terms' that depend on partial derivatives of the metric tensor g ij (x, _ x) with respect to both variables. The 5 If no parametrized path x x(t) is involved, then _ x indicates an arbitrary tangent vector argument. 6 We use Einstein summation convention throughout: Whenever an index occurs twice, once as a lower and once as an upper index, it represents a dummy summation index. Thus e.g.
The reason for considering relative acceleration is that this captures the essential deviation from a constant relative separation velocity observed in the case of Euclidean space. tensor field K i jkℓ (x, _ x) is one of the curvature tensors that can be defined on a Finsler manifold. In the Riemannian case Eq. 30 simplifies considerably, with K i jkℓ (x, _ x) reducing to the so-called Riemann curvature tensor R i jkℓ (x) (a function of x only). In the particular case of interest, with n 3 spatial dimensions, it may be reduced further to involve only the symmetric Ricci tensor field R ij (x), with six local degrees of freedom The interested reader is referred to the literature for technical details (Bao et al., 2000;Shen and Shen, 2016;Rund, 1959).

Riemann-DTI Paradigm: Geodesic Tractography
DTI tensors are represented by symmetric positive definite matrices obtained by fitting their six degrees of freedom at a fiducial point x to the DWI local signal attenuation factor according to the anisotropic Stejskal-Tanner equation (Stejskal and Tanner (1965); Stejskal (1965)), in which q denotes the diffusion sensitizing magnetic gradient direction 8 and b a dimensionful constant related to diffusion time and gradient magnitude used in the scanning protocol. A Riemannian metric may then be constructed proportional to the inverse of the DTI tensor (O'Donnell et al., 2002;Lenglet et al., 2004; The metric tensor is now only dependent on the position x and not on the orientation _ x, so that the Riemannian length of a parameterized curve x x(t) is given by FIGURE 3 | Riemann-DTI paradigm. At a typical point x inside the corpus callosum (black dot in the ROI in the MRI slice) there is a neat alignment of axons (cf. the microscopy image), inducing a diffusion tensor stretched along the preferred orientation as depicted on the right. The ellipsoid is a graphical representation of the quadratic form induced by the Riemannian metric, showing the unit level set g ij (x)y i y j 1 for fixed x, with g ij (x) identified with the (i, j)-th entry of the inverse diffusion tensor. The red arrows illustrate two local tangent vectors. The horizontal one aligns well with the preferred orientation, rendering its Riemannian length relatively short (upper right) as compared to the vertical one of the same Euclidean length (lower right). Geometrically, this length can be inferred via the stack of parallel planes, illustrating the corresponding dual covectors constructed on the basis of the metric via the tangent cone drawn from the tip of each vector to the ellipsoid. In this case one reads off a vector's squared Riemannian length by counting the number of intervals spanned by it: 6, respectively nine for the horizontal and vertical arrows. Note that from a Riemann-intrinsic point of view, i.e. without reference to Euclidean concepts (such as the rigid rotation relating the two arrows), no preferred orientation can be determined. After all, the ellipsoid is, by definition, the Riemannian unit sphere, intended to mask the anisotropic nature of the underlying tissue altogether! 8 The factor b may be absorbed in the tensor q i q j , yielding the so-called B-tensor.
recall Eqs 23, 28. Recall that geodesic tractography entails finding curves x x(t) for which S(x) is locally minimal. This endeavour may be tackled in two different ways, either by directly minimizing Eq. 34, or by solving its Euler-Lagrange equations, Eq. 29. In the Riemannian case at hand, and using an affine ("constant speed") parametrization, these equations reduce to in which the so-called Christoffel symbols of the second kind (Spivak, 1975) are given by linear combinations of derivatives of g ij : Equation 35 may be reformulated in purely geometrical terms by defining a so-called covariant differential operator D (Spivak, 1975), given in terms of the standard differential operator d and Γ-correction terms (a special case of the δ-operator used in Eq. 30), accounting for the non-constant nature of g ij (x). For a vector function v v (t) defined along the curve x x (t) we have so that Eq. 35 is equivalent to For disambiguation of its solution, Eq. 35 requires a pair of initial or boundary conditions The existence of a geodesic constrained by either Eq. 39a or 39b is guaranteed by the Hopf-Rinow theorem (Jost, 2011). For this reason the x-manifold is called geodesically complete. Uniqueness is not provided for the boundary problem, as there may be multiple curves (locally) minimizing Eq. 34 for sufficiently distant endpoints x 0 and x T allowing for multiple potential fiber tracts connecting these endpoints.

Riemann-DTI Paradigm: Geodesic Deviation
Let us consider a single geodesic x x (t) obtained either by minimization of Eq. 34 or by integration of Eq. 35. In order to understand the stability of this geodesic with respect to small perturbations of its initial or boundary conditions, Eq. 39, we may employ perturbation theory in the form of the linear geodesic deviation equations, recall Eq. 30, restricted to the Riemannian case of interest. This allows us to predict the behaviour of neighbouring geodesics without explicitly solving their nonlinear geodesic equations, as long as these are sufficiently close to the fiducial geodesic x x (t), as described in Sengers et al. (2021), in which families of neighbouring geodesics are defined via continuous perturbations of the side conditions, Eq. 39. Here we extend our uncertainty quantification by including the effect of DTI data noise, i.e. perturbations of the metric field g ij (x) along (a full neighbourhood of) a fiducial geodesic, which requires a generalization of Eq. 30 by accounting for an inhomogeneous "force" term. We consider general perturbations of the metric field, as to capture the intrinsic data-induced uncertainty, arising when using the Riemann-DTI tractography paradigm. In practice, these perturbations include perturbations in the DTI tensor and the underlying DWI data.
Thus we consider two types of perturbations, those affecting the intitial or boundary conditions in Eq. 39 and those affecting the metric field, viz. x respectively The formal parameter 0 ≤ ε≪ 1 is dimensionless. We are interested in O(ε) effects; O(ε 2 ) terms are considered irrelevant and will be suppressed 9 .
For an arbitrary geodesic x x (t) satisfying Eq. 35, consider the parameterized family of neighbouring tracks, induced by any of the perturbations in Eqs 40a, 40b, 41. The perturbative, vector-valued function y y (t), defined along the geodesic x (t), must be such that the perturbed path represents itself a geodesic with respect to the metric g ij (x; ε) for any sufficiently small value of ε. Thus the perturbed path, Eq. 42, satisfies Eq. 35 up to O(ε), with the replacement g ij → g ij . The vectors y(t) are referred to as deviation vectors, measuring the extent to which x(t; ε) deviates from x(t) x(t; 0). Substitution of (42) into Eq. 35, with g ij → g ij , yields (cf. Eq. 30) where D/dt denotes the covariant derivative operator along x (t), recall Eq. 37, and R i jkℓ (x) is the Riemann curvature tensor given, in terms of the Christoffel symbols of Eq. 36, by Writing the first order expansion of the Christoffel symbols as Γ akin to Eq. 41, the inhomogeneous term in Eq. 43 is given by 9 For "small enough" y 0 , w 0 , y T , h ij in Eqs 40a, 41 we may set ε 1 without loss of generality In order to appreciate the Eqs 43-46, consider the hypothetical case of a homogeneous, noise free DTI tensor image, inducing a constant metric tensor field. In the absence of noise we have no metric perturbations, so that f i 0, while the x-independent nature of the metric also induces a vanishing Riemann tensor, R i jkℓ 0, as well as a trivial covariant derivative, D/dt d/dt. Equation 43 then reduces to € y i 0, implying a linear evolution of the deviation vector along the geodesic x i x i (t). Both the fiducial geodesic as well as all geodesics induced by this type of deviation are then straight paths by virtue of Eqs 35-38 and Eq. 42. This is what one would expect of a deviation vector in a flat space. In a curved space, with a nonzero Riemann tensor R i jkℓ in Eq. 43, neighbouring geodesics may be locally attracted towards or repelled away from the fiducial one, depending on the nature of curvature. Random data noise, producing a nonzero force term f i in Eq. 43, adds further jitter to the geodesic paths.
Defining the block matrix M ∈ R 6×6 and the column vectors F, Y ∈ R 6 by Equation 43 can be written as an inhomogeneous linear first order system of matrix differential equations This equation may be solved with non-commutative multiplicative calculus (Bashirov et al., 2008;Florack and van Assen, 2012;Gill and Johansen, 1990) using the product integral defined in Eq. 20: The block structure in M, F and Y induces a similar structure in the product integrals: , and (50a) with the help of which the solution y (t) to the initial value problem may be expressed compactly as with Θ(t, σ) Π 11 (t)Ξ 12 (σ) + Π 12 (t)Ξ 22 (σ).
If one considers instead the boundary value problem, the initial tangent vector _ y(0) w 0 must be treated as an unknown and upon setting t T can be solved for in terms of y 0 , y T and Θ(T, σ) f(σ) dσ. The solution of the boundary value problem presents itself as

EXPERIMENTAL RESULTS
So far the theory has been formulated in a continuous setting. In clinical applications of tractography we have a discretized DTI and ditto metric field. This raises the question of sub-voxel tensor field interpolation and discrete differentiation for evaluating Christoffel symbols, recall Eq. 36. These issues are simultaneously handled by using a multiscale representation of the field, Eq. 21, combined with Log-Euclidean interpolation (Arsigny et al. (2006)). For technical details, cf. . Supplied with initial conditions (39a), we may solve for a geodesic by time-integrating Eq. 35 over a specified interval t ∈ (0, T), e.g. using the Runge-Kutta fourth order integration scheme. If instead boundary conditions (39b) are supplied, both direct minimization of Eq. 34 and a nonlinear solver for Eq. 35 are viable solution strategies. Finally, we may discretize the explicit analytical solutions to obtain geodesic deviation vectors from Eq. 53.
To illustrate our perturbative approach, we perform experiments on a DWI dataset from the Human Connectome Project. In our experiments the metric tensor of our Riemann-DTI paradigm is the adjugate of the DTI matrix D, g ij detD D inv ij , cf. . We restrict ourselves to perturbations of the metric given by Eq. 41, while keeping boundary conditions fixed, Eq. 39b i.e. we only consider the influence of perturbations on the DTI tensor field. The effect of boundary variations of type Eq. 40a have recently been considered elsewhere ). These two types of perturbations are invariably present in any tractography algorithm due to DTI domain and codomain uncertainties.
We simulate DTI field perturbations D ij from D ij by a bootstrapping procedure (Whitcher et al., 2008). Corresponding metric perturbations are then determined as h ij g ij −g ij according to Eq. 41 with ε 1, recall footnote 10. Subsequently, the deviation y(t) along an arbitrary track x (t) can be determined by Eq. 53, with y (0) y(T) 0, resulting in a perturbed track x(t) x(t) + y(t). Recall that, apart from discretization, Eq. 53 is an analytical result, rendering this perturbation procedure quite efficient. This data noise simulation procedure is repeated so as to obtain a large collection of perturbed tracks.
For the sake of comparison we explicitly compute, for each simulated metric perturbation g ij , the geodesicx(t) subject to x(0) x(0) andx(T) x(T). We thus obtain two collections of tracks, viz. a set of (computationally intensive) explicitly computed geodesics {x k (t)}, and a set of (analytical, thus computationally efficient) approximations { x k (t)} to geodesics {x k (t)}, where k 1, . . . , N pertains to the track associated with the kth perturbed DTI field. In our experiments we take N 100. Figure 4 illustrates the behaviour of analytical perturbations vs. simulated geodesics. In the first case the analytically computed tracks almost completely coincide with the simulation results, while in the second case a small degree of misalignment is clearly observable. Even though the DTI perturbations are of similar order everywhere, the linear approximation of Eq. 42 apparently does not suffice in all cases, from which we may infer that the mathematical notion of a "sufficiently small" perturbation is case dependent.
In Figure 5 we extend the experiment of Figure 4 to a bundle of geodesics. For each of 4,237 (unperturbed) geodesics starting in the brain stem and ending in the precentral gyrus, their (simulated) perturbed geodesics {x k } and (analytical) approximations { x k } are determined, again with k 1, . . . , N 100, producing two sets of 423,700 tracks each. For each of these sets we construct isosurfaces of the densities obtained by counting the number of tracks passing through each voxel. Subsequently, we threshold the density maps and construct two additional ones, so called difference maps. The difference is set to ±1 for voxels in which fibers are present in only one of the two density maps, and to 0 otherwise. A significant overlap between the analytical and simulated isosurfaces can be observed, confirming the feasibility of our approach.

DISCUSSION
We have presented an extension to the Riemann-DTI paradigm in the form of an uncertainty quantification due to DTI data noise. To appreciate the Riemmanian framework an overview of some of the neccesary tools upon which it is built is given as well.
Variational calculus lies at the basis of geodesic tractography as the mathematical machinery allowing to deduce the geodesic equations from the length functional in the form of Euler-Lagrange equations. The solutions of the geodesic equations are the locally "shortest" paths between two fiducial points in the brain, when "distance" is measured in units of mean free path length by the (inverse) of the DTI tensor. Linear perturbation theory is applied to the geodesic equations, taking into account perturbations in the metric field induced by DTI noise, and leading to an inhomogeneous geodesic deviation equation. Together with multiplicative calculus, notably product integrals, this allow us to express the first order behaviour of the perturbations along the tracks in closed-form. As long as the perturbations in the DTI tensor are sufficiently small (which is case dependent), the analytically computed perturbed tracks do not significantly differ from the explicitly computed tracks, which are geodesics with respect to the perturbed DTI metric. This is confirmed by our experiments, both in the case of a single (unperturbed) geodesic as well as for bundles of geodesics.
In the experiments we compared both the analytically computed perturbations as well as the simulated geodesics with respect to the perturbed metric. In practice, however, one may prefer the analytical approach, as it offers a strong computational advantage. By doing brute force simulations we need to solve non-linear differential equations to find the perturbed tracks, while in our analytical approach we employ a linear differential equation with a closed-form solution that can be efficiently computed.
We recall that geodesic completeness ensures (at least) one connecting track for any pair of sufficiently distant endpoints. In case of multiple connections external knowledge may be necessary to select the most plausible one from the biological point of view. For example, tractometric features can be used to classify and subsequently prune connections according to their diffusion-related characteristics (Colby et al. (2012);De Santis et al. (2014)). We will address this important point in future work.
Our perturbative analysis may be extended past the Riemannian paradigm by using Finsler geometry, which appears to be a natural extension for general DWI models beyond DTI. This extension removes the quadratic restriction in the metric and allows for a more accurate modelling of the DWI signal capturing additional features such as intra-voxel heterogeneity, e.g. due to partial volume effects or crossing fibers. The Finslerian approach yields stronger descriptive power, at the cost of a mathematically similar but more The red surface contains is the intersection of the analytical (green) and simulated (purple) isosurfaces. The blue surface indicates the region which contains analytically perturbed tracks but not simulated tracks, and vice versa for the yellow surface. The background DWI image is included purely for reference to the location of the bundle in the brain. However, since the slice of the DWI image is offset from the bundle itself, this visualisation should not be used to judge the quality of the tracks, merely that of the difference between the analytically computed perturbations and their simulated counterparts.
Frontiers in Computer Science | www.frontiersin.org October 2021 | Volume 3 | Article 718131 cumbersome representation. The Finslerian counterpart of the Riemannian geodesic deviation equation may likewise inform us on the effect of perturbations on general DWI data.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: ConnectomeDB repository on https://db. humanconnectome.org; dataset "WU-Minn HCP Data-1200 Subjects"; subject 100307.

AUTHOR CONTRIBUTIONS
RS and LF contributed equally with different emphases. Abstract and Sections 1 and 2 were written by LF. Theory outlined in Section 3 is joint work by RS and LF. Experiments described in Section 3 are conducted by RS, including algorithmic implementations. Section 4 is joint work by all authors. AF contributed to the overall setup, critical proofreading, and editorial, and to the fundamental research on which this article is based.