Normal and equivolumetric coordinate systems for cortical areas

We describe coordinate systems adapted for the space between two surfaces, such as those delineating the highly folded cortex in mammalian brains. These systems are estimated in order to satisfy geometric priors, including streamline normality or equivolumetric conditions on layers. We give a precise mathematical formulation of these problems, and present numerical simulations based on diffeomorphic registration methods, comparing them with recent approaches. Our method involves• Diffeomorphic registration of inner and outer folded folded surfaces.• Followed by equivolumetric reparametrization of layers to yield coordinate system.


Introduction
Anatomical regions of interest extracted from 3D biomedical imaging data often appear as volumes separated by "upper" and "lower" surfaces.These include, in particular, the highly folded cortical regions or areas of the mammalian brain.This paper focuses on methods that parametrize such volumes, using coordinate systems that are naturally aligned with the encompassing surfaces.Our contributions are, on one hand, to formalize a notion of laminar coordinate systems and, on the other hand, discuss within this formalism the concept of equivolumetric coordinates [5,6], providing an interpretation of two similar methods [21,16] and introducing a new one based on diffeomorphic registration methods.
The elegant laminar structure of a typical cortical area [8,22] is summarized as follows.The folding of the area serves to maximize its surface area in a confined cranial space.The neural tissue (grey matter) within the area contains mostly neuronal cell bodies and unmyelinated fibers.Cortical areas (which number in the hundreds in the human brain) are connected via white matter containing axonal, usually myelinated, fibers.Each cortical area is composed of Figure 1: Idealisation of Bok's equivolumetric hypothesis for the layers in the folded cortex showing that the deep layers (V-VI) are thin in the sulcus and thick in the gyrus respectively characterized by regions of positive and negative inward curvature and conversely for the upper layers (I-III) fundamental units called cortical columns that traverse vertically from the white matter to the surface just below the pial matter.Finally, the cortical area is composed of six layers which are stacked horizontally on top of each other.Bok [5,6] observed that to maintain the laminar structure in highly folded regions that thin layers in one part became thicker in another part.This observation led to the hypothesis that the cortex satisfies an equivolumetric property, illustrated in Figure 1, which provides the theoretical motivation of this paper.
The rest of the paper is as follows.Section 2 introduces formal definitions of laminar coordinates and describes methods for constructing them.Section 3 describes how Bok's equivolumetric hypothesis can be implemented.Section 4 presents numerical simulations.
2 Laminar Coordinate Systems and Thickness

Notation
We introduce some mathematical notation in order to describe 3D coordinate systems parametrizing an open space between two surfaces.We will call "laminar coordinate system" (see Figure 2 for an example) a special case of a foliation of this open set, with two special leaves provided by the two surfaces.
More precisely, let S 0 and S 1 be 2D submanifolds (surfaces) in R 3 , such that Figure 2: A laminar coordinate system between two surfaces, represented by the red scaffold.
The surfaces S t = ψ(t, S 0 ) are the leaves of the foliation, and will be referred to as "layers" while the curves γ x (t) = ψ(t, x), t ∈ [0, 1] will be referred to as "streamlines."These layers and streamlines are shown as blue surfaces and vertical red lines respectively in Figure 2. The set Ω = ψ([0, 1] × S 0 ) is the closure of an open subset Ω of R 3 defining the space between the two surfaces.
In addition, a coordinate system defined as such provides a definition of "thickness" of Ω at x ∈ S 0 , given by the length of the streamline starting at x, namely,

Coordinate Systems
There are several methods for building coordinate systems.Three approaches are described as follows.

Level Set Methods
Assume that Ω ⊂ R 3 and that its boundary has two connected components, providing S 0 and S 1 , which must therefore be closed surfaces.The typical application is when Ω is the region between two non-intersecting nested spherelike surfaces.
Assume that one is given a C 2 submersion F : Ω → [0, 1] such that S 0 = F −1 (0) and S 1 = F −1 (1).(This assumes that F is onto and that ∇F (x) = 0 for all x ∈ Ω.) Define the vector field and let ψ(t, x), t ∈ [0, 1], x ∈ S 0 satisfy with ψ(0, x) = x, so that ψ is the restriction to S 0 of the flow associated with the vector field v.Then, one has which implies that F (ψ(t, x)) = t for all x ∈ S 0 and t ∈ [0, 1].The function ψ : [0, 1] × S 0 → Ω provides a laminar coordinate system of Ω as defined above.The thickness along streamlines is in particular defined by Notice that, since the layers S t coincide with the level sets of F , the vector field v (and therefore the streamlines) are necessarily perpendicular to them.
Given S 0 and S 1 , one therefore needs to define a suitable function F .Jones et al. [12] defined F as the solution of the Laplace equation ∆F = 0 with boundary conditions F = 0 on S 0 and F = 1 on S 1 .Similarly, Waehnert et al. [21] proposed a construction starting with level-set representations of S 0 and S 1 , represented by two functions U 0 and U 1 such that S i = U −1 i (0), i = 0, 1.This leads to the definition where S is a smoothing operator based on a topology-preserving mean-curvature motion equation.More precisely, given a function U targ , the function SU targ is obtained as the limit when time s → ∞ of the solution U (s, •) of

Diffeomorphic Volume Mapping
Das et al. [9] assumed that the surfaces S 0 and S 1 are represented as boundaries of two open set Ω 0 and Ω 1 (so that S 0 = ∂Ω 0 and S 1 = ∂Ω 1 ) with Ω0 ⊂ Ω 1 .Then a variant [2] of the large deformation diffeomorphic metric mapping (LDDMM) algorithm [3] was used to estimate a flow of diffeomorphisms ϕ(t, •) : The function ψ can then be defined as the restriction of ϕ to [0, 1] × S 0 .

Surface Mapping with Normality Constraints
The foregoing methods are volumetric: they work on 3D regions rather than directly on surfaces, which are represented either as level sets or as boundaries of open sets.They are, in particular, not well adapted to analyze regions delimited by open surfaces.So following [18], a version of LDDMM adapted to surface mapping was combined with normality constraints in order to estimate laminar coordinate systems.This is now described in some detail.
Assume that S 0 is parametrized over a bounded open subset M ⊂ R 2 (or more generally over a 2D manifold M with or without boundary) in the form S 0 = q 0 (M ), where q 0 is an embedding of M into R 3 .The LDDMM surface registration algorithm solves an optimal control problem minimizing, over all time-dependent vector fields in a reproducing kernel Hilbert space V , subject to q(0) = q 0 and ∂ t q(t) = v(t) • q(t).Here, D is a reparametrizationinvariant discrepancy measure between (unparametrized) surfaces.Several versions of this cost function have been introduced, based on representation of surfaces as currents [20], varifolds [7,13] or normal cycles [19].Assume that S 0 and S 1 are triangulated surfaces and that the cost function D is replaced by a discrete approximation, still denoted D.Then, the optimization problem can be reduced to one tracking explicitly the evolution of the vertices of the triangulation, using the reproducing kernel of V denoted as K.This kernel is a matrix-valued function of two variables x, y ∈ R 3 such that, for all α, y ∈ R 3 , the vector field x → K(x, y)α belongs to V and for all v ∈ V , where the left-hand side denotes the inner product in V .
Denote as q 0 = (q 0 (1), . . ., q 0 (N )) the vertices of S 0 .The reduced problem is expressed in terms of evolving vertices q(t) = (q(t, 1), . . ., q(t, N )) and vectors α(t) = (α(t, 1), . . ., α(t, N )), t ∈ [0, 1], minimizing (letting S(t) denote the triangulated surface with vertices q(t) and same topology (faces) as S 0 ) for k = 1, . . ., N .Moreover, the optimal vector field at time t is given by ( Figure 3: Synthetic data for thickness estimation between an inner ring of fixed radius and an outer one of variable radius.The colors on the outer ring show thickness values increasing from blue to red. The interest of this formulation is that the trajectories t → q(t, k) for k = 1, . . ., N directly provide the streamlines starting from the vertices q 0 (k), k = 1, . . ., N of the triangulation of S 0 .This algorithm is modified by imposing a constraint ensuring that these streamlines are perpendicular to the evolving layers [18].In the continuous setting, where, for each t ∈ [0, 1] q(t, •) is defined on the manifold M , this constraint can be formulated as v(t, q(t, s)) = λ(t, s)ν S(t) (q(t, s)), t ∈ [0, 1], s ∈ M .Here ν S (x) denotes the (positively oriented) unit normal to an oriented surface S at x ∈ M , and λ(t, s) is a scalar that we assume to be non-negative (with a proper orientation of S(t) to prevent the trajectories from backtracking).The constraint is implemented in the equivalent form that is discretized on the (evolving) triangulation of S(t).The resulting constrained optimization problem is solved using an augmented Lagrangian method [17] with each gradient descent step implemented using a limited-memory BFGS method (for details of methods using LDDMM with constraints, see [1]). Figure 3 illustrates a synthetic example of how the thickness map can be estimated using this method.
Note that the LDDMM algorithm implicitly provides a flow of diffeomorphisms on the whole space R 3 , given by solution of Even if we will not use it in the following, it is interesting to notice is that this method generally provides a level-set formulation of the laminar coordinates.Indeed, under the mild assumption that v(t, q(t, s)) never vanishes, the function ψ : [0, 1] × S 0 → R 3 defined by ψ(t, q 0 (s)) = q(t, s) is an immersion.In practice, this mapping is in addition one-to-one and one can define without ambiguity a scalar function F on a domain sandwiched by the surfaces S 0 and S 1 (or, more precisely, by S 0 and S(1) S 1 ),by for s ∈ M .By construction, the streamlines are perpendicular to the level set of this function.

Strict Condition
As mentioned above, Bok's hypothesis requires that the cortical layers satisfy an equivolumetric constraint.Define a "cortical tube" as the volume delimited by the cortical columns stemming from the inner (grey-white) matter surface to the outer (pial) surface.The hypothesis requires that the volume delimited by the intersection of cortical tubes and cortical surfaces remains roughly constant when the base patch is "translated" along the inner surface.We want to formalize this into a definition of "equivolumetric laminar coordinates." Consider a laminar coordinate system ψ : [0, 1] × S 0 → R 3 .For x ∈ S 0 , consider a small surface element δS 0 located at x and the infinitesimal tube ψ([0, 1] × S 0 ).Introduce a local chart m : U → δS 0 on δS 0 , where U is an open subset of R 2 .Let ψ m (t, α, β) = ψ(t, m(α, β)).Then the volume of the tube between layers t 0 and t 1 is given by Here, we have used the notation ν(t, x) = ν S(t) (ψ(t, x)).We also have denoted by dvol S0 the volume form on S 0 , which is given, in the local chart, by |∂ α m × ∂ β m|dαdβ, and by σ(t, x) the surface Jacobian induced by ψ(t, •) : S 0 → S(t) (infinitesimal ratio of area), defined in the chart by One has the following result that describes the evolution of σ as a function of t.  with ρ(y) ⊥ N (y) for all y ∈ Ω.Let ρ t denote the restriction of ρ to S(t).One has where div S(t) is the divergence operator on S(t) and H S(t) the mean curvature on the same surface.
The notation used in this proposition is illustrated in Fig. 4. Recall that the divergence operator of a vector field ρ on a surface S is (for p ∈ S) where e 1 , e 2 is any orthonormal basis of T p S (the tangent space to S at p).
Using the same basis, the mean curvature is given by When convenient, will use the notation: H(t, x) = H S(t) (ψ(t, x)) to represent the mean curvature along a streamline.
Proof.We make the computation in a chart m : U ⊂ R 2 → S 0 and let One has ).These identities can be proved from ( 4) and ( 5) by introducing an orthonormal basis (e 1 , e 2 ) of T ψm(t,x) S(t) and expanding the left-hand sides as functions of the coordinates of e α and e β in this basis (cf.[24], lemma 3.19).Using this, we get from which one deduces that Define the "equivolumetric thickness" along the streamline starting at x by A strict interpretation of Bok's hypothesis requires that this expression does not depend on x, i.e., that there exists a function t → λ(t) such that for all x ∈ S 0 .(In which case γ(t, x) = t 0 λ(u)du.)One can, without loss of generality, assume that λ does not depend on t, which can be achieved by applying a time change to the evolution.More precisely, one can replace ψ by ψ such that ψ(τ (t), x) = ψ(t, x) with Therefore assuming that λ is constant, we now apply Proposition 1 to obtain a surface propagation equation that is equivalent to (7).We will decompose as required by (7).One then has the following proposition.
Proposition 2. A laminar coordinate system ψ satisfies Bok's hypothesis if there exists a vector field ρ tangent to the layers defined by ψ and a constant λ with ψ(0, x) = x and σ(0, x) = 1 for all x ∈ S 0 .
Equation ( 8) provides an evolution equation controlled by ρ (such that at all times, ρ(t, •) is a vector field on the evolving surface S(t)) whose solution satisfies the equivolumetric hypothesis.A detailed study of this system of equations (including its well-posedness for a given choice of ρ), and of the optimal control problem consisting in optimizing ρ with the constraint that S(1) = S 1 are challenging open problems that will not be addressed in this paper.Even if possible, it would also be counter-intuitive to require a constant equivolumetric thickness, γ(1, x), in (6).So in the next section, we discuss a solution to a simpler problem that we refer to as a "localised" Bok's hypothesis.

Localised Bok's hypothesis
In this section, we replace the strong constraint ( 7) by a weaker one for a given function c 0 .Here again, there is no loss of generality in assuming that λ is constant, and, including if needed this constant in c 0 , in taking λ = 1 (so that c 0 coincides with the equivolumetric thickness γ(1, x)).Proposition 2 can be directly extended in this setting.
Proposition 3. A laminar coordinate system ψ satisfies the localised form of Bok's hypothesis for a given equivolumetric thickness c 0 if there exists a vector field ρ tangent to the layers defined by ψ such that with ψ(0, x) = x and σ(0, x) = 1 for all x ∈ S 0 .
Because the function c 0 can be chosen freely, (9) (or a solution of system (10)) is much easier to obtain while satisfying the condition ψ(1, S 0 ) = S 1 and we now show that it can be achieved starting from any laminar coordinate system ψ by applying a space-dependent time change.Indeed, let ψ(τ (t, x), x) = ψ(t, x), where, for each x ∈ S 0 , t → τ (t, x) is an increasing differentiable function from [0, 1] onto [0, 1].Then, introducing as above local coordinates α and β on S 0 , we have As a consequence: In order that (9) holds for ψ, we therefore need to define τ so that (for some function c 0 ) In order to have τ (0, m) = 0, τ (1, m) = 1, we need and then So, the time change is provided by the relative volumetric depth the streamlines (that are left unchanged in the operation).The equivolumetric layers at level are provided by points ψ(t, x) along the streamlines satisfying τ (t, x) = (in the original parametrization).

Interpretion of recent models
We now interpret two recent attempts to model Bok's hypothesis [21,16] in our framework.
Waehnert et al. [21] start with streamlines estimated using (1).Equivolumetric layers are then estimated using the same equation, replacing the constant value ρ by a function ρ(x), that is determined as follows.Using our notation, one first makes the assumption that therefore making a linear approximation of the surface change.Here, one takes t = s/θ, where s is the arc length along the streamline and θ is the thickness (the length of the streamline).For x ∈ S 0 , the value of σ(1, x) is estimated as a function of the curvatures at both ends of the streamline starting at x (we refer to [21,14] for more details and justification).Integrating along streamlines, which are perpendicular to the layers because of the level set formulation, one obtains an expression of the equivolumetric depth given by Given ∈ [0, 1], one defines a target level ρ (x) corresponding to the layer at equivolume by solving V x (ρ) = αV x (1), which is a quadratic equation in ρ.

Numerical Implementation
Because they rely on level sets, the two recent approaches [21,16] are Eulerian, i.e., they work in the 3D volume Ω and layers are isosurfaces associated with scalar functions defined on Ω while streamlines are integrated using the gradient of these functions.
Our approach is different in that it is directly modeling layers as parametrized surfaces S(t) = ψ(t, S 0 ), so that our implementation is based on a triangulation of S 0 and a discretization of the time interval.The streamlines are obtained as solutions of the ODE ∂ t y = v(t, y), where v is obtained from the LDDMM algorithm and provided by (2).Importantly, this time-dependent vector field is discretized in time only, and known analytically as a function of y.In particular, its space derivatives can be evaluated without approximation.It also specifies a flow of diffeomorphisms of R 3 through the equation with ϕ(0, x) = x for all x ∈ R 3 .
The surface Jacobian σ can be evaluated using the evolving triangulated surfaces: if x is a vertex on S 0 , we let a(0, x) denote the area of the one-ring centered at x (the union of all triangles that contain x).Similarly, we let a(t, x) denote the area of the one-ring around ϕ(t, x) in the triangulated surface S(t) = ϕ(T, S 0 ) (which has the same triangle structure as S 0 ).On can then define σ(t, x) = a(t, x) a(0, x) .
This is the approximation that are used in our simulations, and it is accurate provided that the triangulation of S 0 is fine enough without flat triangles.An alternative procedure is also possible, since one has, in this context for x ∈ S 0 .(The "−T " exponent refers to the inverse of the transpose matrix.)Recall that ϕ(t, •) is (for fixed time t) a diffeomorphism of R 3 , hence defined on the whole space (unlike ψ, which, for laminar coordinates, is only defined on S 0 , and in this special case, is defined as the restriction of ϕ to this surface).This implies that ∂ x ϕ(t, x) is a 3 × 3 matrix.To prove (12) one can just notice that, and use the fact that for any matrix A and vector u and v, one has The time evolution of the vector ζ(t, x) = det(∂ x ϕ(t, x))∂ x ϕ(t, x) −T ν(0, x) is provided by This can be integrated along streamlines using the expression of v in (2), from which, as mentioned, space derivatives can be evaluated exactly.

Results
Results of the numerical implementation are presented for three cases.Figure 5 shows the result for the synthetic data of Figure 3. Figures 6 and 7 respectively show the results for the marmoset auditory cortex (obtained from [23]) and feline auditory cortical regions (obtained from [4]).Here, the proposed method is compared with those of Waehnert et al. [21] and Leprince et al. [16] computed via Github packages [15,11].Figure 8 shows the corresponding cumulative distribution of distances of equivolumetric surfaces at t = 0.25, 0.5, 1.0 relative to those via the proposed method.Here, for surfaces S 1 and S 0 with vertices x i ∈ S 1 and y j ∈ S 0 , the distance at the i th vertex of S 1 is where n = arg min j d(x i , y j ) and m = arg min k d(x k , y n ).This FreeSurfer distance [10] returns a value for every vertex and making it more robust against outliers that arise with the Hausdorff distance.
Figure 5: Estimated equivolumetric layers (in red) at times t = 0.3 and 0.6 for the synthetic example of Figure 3 .

Discussion
A unified theoretical framework for describing laminar coordinate systems for cortical regions has been developed.Herein, several algorithms including an interpretation of volumeteric-based ones are offered such that equivolumetric layers consistent with Bok's hypothesis can be computed.Granted that in the marmoset the auditory cortex is the only area that is folded with a gyral crown, the protrusions of the upper layers seen with the method of Leprince et al. [16] may be attributed to the diverging normal vector fields as one approaches the outer surface.This divergence warrants finer discretization which may be computationally expensive.The primary and higherorder auditory cortical regions in the cat reveals greater differences between the three methods.The t = 0.25 layer appears to be closer to the sulcal fundi for Leprince et al. and our methods; in contrast the t = 0.75 layer appears to be closer to the gyral crowns for Waehnert et al. and our methods.These differences can be quantified via a distance metric (Figure 8) and may be attributed to the representation of curvature -direct or indirect-in the computations.Future work will examine the how disease and disorder affect the equivolumetric depths of layers, pyramidal cells and other cortical elements in 3D.This will build upon previous work in 2D [4].Our laminar coordinate system also has a straightforward application to cortical layer segmentation.As a prior, it could Figure 6: Equivolumetric layers for a marmoset auditory cortex at t = 0.25 (orange) and t = 0.75 (yellow) using Laplacian [16,15] (top), level set [21,11] (middle) and the proposed method (bottom).
The time change introduced to reparametrize the coordinates in order to make them compliant with Bok's hypothesis left the streamlines invariant while changing the layers.As a consequence, if streamlines were perpendicular to the layers to start with, this property is generally lost after reparametrization.Finding an equivolumetric coordinate system with perpendicular streamlines is a significantly more arduous problem.Using Proposition 3, in which one must set ρ = 0, one sees that this problem requires to estimate a scalar field c 0 on S 0 such that the solution of (10) satisfies ψ(1, S 0 ) = S 1 .Whether this inverse problem is well posed, and whether stable numerical algorithms can be designed to solve it, are open questions that we plan to address in future work.

Figure 4 :
Figure 4: Illustration of the notation used in Proposition 1.The evolution of intermediate surfaces (red) is captured by a vector field w that decomposes into a tangential part ρ and a normal part ζN .