Exact analytical PGSE signal for diffusion confined to a cylindrical surface using a spectral Laplacian formalism

1 Exact analytical PGSE signal for diffusion confined to a cylindrical surface using a spectral Laplacian formalism Erick J Canales-R odríguez 1,2,3,* , Chantal M.W. Tax 4,5 , Juan Manuel Gó rriz 1,2,3 , Derek K. Jones 5 , Jean-Philippe Thiran 6,7,8 , Jonathan Rafael-Patiño 6,7 1 Department of Signal Theory, Networ king and Commun ications, Uni versity of Granad a, Granada, Spain 2 Andalusian Research Insti tute in Data Science and Co mputational Intelligenc e (DaSCI), Uni versity of Granada, Spain 3 Research Centre for Inf ormation and C ommunication Technologies (CITIC) , University of Granada, Granada, Spain 4 Image Sciences Insti tute, Universit y Medical Cent er Utrecht, The N etherlands. 5 Cardiff University Brain R esearch I maging Centre (CU BRIC), Cardiff Universi ty, Cardiff, Wales , United Kingdom. 6 Signal Processing Labora tory 5 (LTS5), Ecole Polytechniq ue Fédérale de Lausanne (EPFL), Lausanne, Switzerland. 7 Computational Medic al Imaging & Machine Learn ing Section, Cent er for Biomedical Imaging (CIB M), Lausanne, Switzerlan d 8 Department of Radiol ogy, Centre Hospitali er Universi taire Vaudois (CHUV) , Lausanne, Switzerlan d * Corresponding auth or: ejcanalesr@ug r.es 2 Abstract Pulsed- gr adient spin-echo (PGSE) MRI experiments probe molecular self-diffusion thr ough spin phase accumulation under time-depend ent magnetic field gradients. For diffusion confined to cylindrical surfaces, existing analytical sign al models typically rely on the narrow-pulse limit, approximate treat ments of finite gradient durati ons, or the Gaussian phase ap proximation, which bec ome increasingly inaccurat e at hig h d iffusion weightings. Here, we derive an exact analytical s olution of the Bloch – Torrey equ ation for diffusion confined to a cylindrical surface under finite PGSE gradien ts and obtain the corresponding diffusion MRI signal expression valid for arbitrary gradien t durations and separations. The derivation is based on a spectral matrix f ormalism of th e Laplace op erator in the eig enbasis o f the confinin g g eometry. The signal is expressed as a product of non-commutin g matrix exponentials, without approximations to the diffusion propagator or the spin phase distribution. We furth er introduce a reduced real spectral basis exploiting the symmetry of the cylindrical surface, sub stantially improving computational efficiency. Building on this exact f ormulation, we de velop efficient numerical strategi es for rep eated signal evaluations, including a Strang splitting approximat ion of the m a trix exponentials and an efficient computation of the spherical mean signal using Gauss – Legendre quadrature. The analytical signal is validated against Mo n te Carlo simulations o v er a wide range of cylinder radii and experimental parameters. The accelerated implementations are benchmarked against the exact formulation to quantify accuracy – runtime trade-offs. These results est ab lish a computationally efficient framework for evaluating directional and orientationally averaged diffusion M RI sign als in applications req uiring large n umbers of model evaluati ons. Keywords : Diffusion MRI; Laplace operator; Pulsed‐grad ient spin‐echo sequen ce; Cylindrical surface; Monte Carlo diffusi on simulations ; Myelin sheath radi us. 3 1. Introduction Magnetic Resonance Imagi ng (MRI) e xperiments sensitized to molecular self-diffusion provide a powerful framework to probe transport processes in confined environments and to infer geometric features of restricting boundaries in porous media and biologica l tissues 1 – 4 . In pulsed-gradient spin-echo (PGSE) experiments 5 , diffusion sensitization arises from th e accumulation of spi n phas e under time-dependent magnetic field gradi ents, such that the measured macroscopic signal reflects ensemble-averaged properties of the underlyin g stochastic motion. Establishing rigorous links between the measur ed signa l and the geometr y of confinement re mains a central pr o b lem in the the ory of diff usion MRI (dMRI) 6 – 18 . From a the oretical perspe ctive, the dMRI signal from water molecules in a bounded domain can be described exactly by solving the Bl o ch – Torrey equation. Spectral and matrix-based approaches employing the eigenfunctions of the Laplace o pera to r have been developed to address this problem 19 – 28 . These frameworks provide a rigorous description of restricted diffusion and, in principle, allow one to treat arbitrary gradient wav eforms and boundary condition s 29,30 . However, for specific geo metries of interest, practical signal models often adopt simplifying assumptions (most notably the narrow -p ulse limit, approximate treatments of finite gradient durations 31 , or the Gaussian-phase approximation 32 – 34 ) in order to obtain expres sions that ar e analytically trac table and co mputationall y efficient. Diffusion confin ed to cylin drical surfaces constitutes a particularly relevant cas e. In biol ogical tissues , this geometry has been proposed as an effective model for myelin water diffusion, in which water molecules are c o nstrain ed to move along the surface of the myelin sheath surround ing axons 35,36 . This assumption is support ed by the ultrastructure o f myelin, where adjacent lamellae ar e separ ated by narrow aque ous gaps with an effective thickness of approxim ately 3 – 4 nm, as inferred from diffraction- and scattering- based reconstructions 37,38 . In a previous study, we derived an exact analytical expression for the PGSE signal arising from diffusion confined to a cylindrical surface in the narrow-pulse lim it and proposed an approximate extension to account for finite gradient durations 35 . While that formulation showed good agreement with M ont e Carlo simulations over a broad range of experimenta l conditions, noticeable discrepancies emerged for large diffusivities and relatively high diffusion weig htings, highlighting the limitations of the ap proximations when fin ite gradient pulse effects bec ome signi ficant. The main objective of the p rese n t work is t o res olve th is limitation by d eriving an ex act analytical solution of the Bloch – Torrey equation for diffusion taking place on the two-dimensional cylindrical manifold , without invoking any approxi mation o n the diffusion propagator, resulting spin phase distribution, or gradient pulse duration. Th e derivation is based on the spectral matrix f ormalism of th e Laplace operator 29,30 , specifically constructed for the cylindrical surface geometry. By explicitly accounting for the non- commutativity of the o per ators involved during the different segments of the PGSE sequence, we obtain a clo s ed-fo r m dMRI signal expression valid for arbitrary combinati ons of rectangular gradient duration and separation. The paper is organized as follows. Section 2 outlines the theoretical framework for describing confined diffusion in PGSE experi m e nts using the Laplacian spectral matrix f ormalism. In Section 3, this fra mework is applied to diffusion confined to a cylindrical surface, and the resulting analytical signal expression is derived. Building o n this exact formulation, we also develop efficient numerical strategies for repeated signal evaluations, including a Strang splitting approximation and an efficient computation of the spherica l mean signal using Gauss – Legendre quadrature. Details of the Monte Carlo diffusion simulations used for numerical validation ar e also pr ovided. Secti o n 4 presents the comparison betwee n analytical and Monte 4 Carlo signals and bench marks the accelerated i mplementati ons against the exact formulati o n to quantify accuracy – runtime trade-off s. The scope and limitations of the proposed approach are discussed in Section 5. Additional technical d et a ils o f the mathema tical derivations are provided in th e Appendices. 2 Theory In this section, we briefly recall the theoretical framework for describin g confined diffusion based on the spectral decomposition of the Laplace operator. The use of Laplace operator eig enfunctions to solve the Bloch – Torrey equation was originally introdu ced b y Robertson 19 , and subsequently e xtended throu gh a variety of analytical and numerical approa ches. N otable de velopments include th e multip le narr ow-pulse approximation 20 , matrix-based formulati o n s enabli ng efficient numerical i mplementations 21,22 , and additional spectral treatments 24,39 . A more complete theoretical description was later formulated by Axelrod and Sen 25 and further generalized and reformulated within the multiple co rr elation function framework by Gr ebenkov 28 . I n what follows, we summarize the key relations o f this spectral matrix approach that are required for the pre sent work, before applying it to diffusion confined to a cylin drical surface. 2.1 Bloch – Torrey equ ation an d spectral approa ch The Bloch-Torrey equation in the rotating frame of reference, d escribing the time ev olution of the complex-valued transverse magnetization ( ) , Mt r , ari sing from spins diffusing with diffusion coefficient D in a confined/b ounded domain  is 40 ( ) ( ) ( ) ( ) ( ) 2 , , ,, Mt D M t i G t M t t    =  −   r r g r r , (1) where ( ) Gt is the magnitude of the applied time-depende nt linear magnetic field grad ient pointing along the unit vector g , r is the spin’s spatial position vector,  is the gyromagne tic ratio of the nuclei, and 2 ,   denotes the Laplace oper ator 2  o n the domain  r where the diffusion process takes place , which depends on the boundary conditi ons describing how th e spins interact with the boundaries (e.g.,  ={ Neumann, Dirichl et, Robi n}). In practice, it is convenient to represent ( ) , Mt r in the orthonormal basis of eigenfunction s   n  resulting from th e spectral dec omposition of the Laplac e opera tor 2 ,   . These complex-valued eigenfunctions are determined by the f ollowing eigenvalue proble m 29 : ( ) ( ) 2 , , n n n     − = rr (2) where dependin g o n the geometry of the problem, n belongs to the set of natural numbers   0 , 1 , n  =  or t o the set of in tegers   , , 1 , 0 , 1 , n  = −  −  . The ass o ciat ed eigen values ar e 5 ordered as 01 0     when n  , or 1 0 1    −     when n  , and the orthonorm al eigenfunctions satisf y the propert y: ( ) ( ) * , n m n m d     =  r r r , (3) where * m  is the compl ex-conjugate of m  and nm  is the Kronecker del ta: it is equ al to 1 if nm = , and 0 if nm  . In this basis, the magnetization can be e xpanded as: ( ) ( ) ( ) ,, nn n M t c t  =  rr (4) where ( ) n ct are time-dependen t fu nctions to be esti mated. Substituting Eq. (4) in to Eq. (1), and usin g the definitio n in Eq. (2), we obtain ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) . n n n n n n n n n n ct D c t i G t c t t       = − −      r r g r r (5) This relationship can be simplified by multiplying both sides of the equati o n by * m  , then integra ting across the spatial variable and usin g the identit y in Eq (3), which results in ( ) ( ) ( ) ( ) ( ) ( ) ( ) * . n n n n n m n ct Dc t i G t c t d t       = − −     r g r r r (6) By defining a column vector of time-dependent functi ons, i.e., ( ) ( ) ( ) 01 ,, T t c t c t =     c when n  , or ( ) ( ) ( ) ( ) 1 0 1 , , , , T t c t c t c t − =     c when n  , we obtain the following system of equations in matrix-form: ( ) ( ) ( ) ( ) , t D t i G t t t   = − −  c Λc Bc (7) where Λ is a diagonal matrix of the eigenvalues estimat ed from Eq. (2) :   ( ) 01 ,, d ia g  = Λ when n  , or   ( ) 1 0 1 , , , , d ia g    − = Λ when n  . The entries of matrix B are determined b y solving the integra l ( ) ( ) ( ) * . nm n m Bd   =  r g r r r (8) The solution to the system of differential equations in Eq. (7) is given by 24,39 6 ( ) ( ) ( ) ( ) 0 0, T D i G t dt te  −+  = ΛB cc (9) where e  denotes the time-ordered matrix exponential an d ( ) 0 c is the vector evaluated at time z ero. In next sections we will determine matrices Λ and B for the confining domain under study, but first, we will relate the dMRI signa l with the time dependent vector ( ) t c and the eigenfunctions in the next subsection. 2. 2 Diffusion MR I signal: Expa nsion on the spectral b asis The measured dMRI signal is obtained by spatiall y int egrating the transverse magnetization, which can be represented in the basis of eigenfu nctions introduced i n the previous section: ( ) ( ) ( ) ( ) ( ) , , nn n E t M t d c t d t    = = =    rr rr wc ( 10 ) where w is the v ector with entries equal to the spatial integral of the eigenfunctions: ( ) ( ) 01 ,, dd    =     w r r r r when n  , or ( ) ( ) ( ) 1 0 1 , , , , d d d    −     =       w r r r r r r when n  . This expression esta blishes a direct rela tionship betw een the dMRI sign al and vector ( ) t c . 2.3 Analytical solut ion for the PGSE experiment In a PGSE sequence wi th tw o rec tangular pulses of du ration  , separa ted b y  − , the gradient function can be represent ed as piecewise constant in 3 segm ents 5 : 1. First pulse: ( )   f or 0, G t G t  = , 2. Free evolution : ( )  ( 0 f o r , G t t  =   , 3. Second pulse: ( )  ( f or , G t G t  = −    + . For this sequence, the general s olution given in Eq. (9) becomes 24,29,39 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 0, 0. D dt D i G t dt D i G t dt D i G D D i G t e e e e e e            +  − − + −+ − − −  − − +   = = Λ Λ B ΛB Λ B Λ Λ B cc c ( 11 ) It should be noted that this relationshi p involves matrix exponenti als applied to square matrices, which are not equal t o the ordina ry exponential functi o n. As the involved matrices do not co mmute, we cannot regroup them in to a single exponential. The n umerical evaluation must be perfor med b y multiplying the 7 exponential matrices from right to left: propagating the effect of the first gradient pu lse on the initial state, then including the effect of the free e volution part on the resulting vector, and finally the effect of the second pulse. Neverthe less, as matri x Λ is diag onal, the matrix e xponential in the second term in Eq. ( 11 ) can be effici ently comp uted by exponentiatin g each entry on the main diag onal. 3 . Methods 3.1 Diffusion con fined on cylind rical surfaces: Compu ti n g matrices Λ and B After having reviewed the general theory o f the spectral approach and resulting dMRI signal, the last missing ingredient of the recipe is to compute the matrices Λ and B for our problem of interest: for diffusion confined to a c ylindrical surface. In this section, we will focus on deriving an ex p ression for the dMRI signal E ⊥ fo r a diffusi on gradient applied perpendicu lar ly to the cylinder’s main a xis, which thus, depend s on the spin’s displace ments on the 2D circumferenc e of the circle perpendicular to th is axis. In polar coordinates, the Laplace op erator for this geometry reduces to its angu lar component part : 2 2 , 22 1 r    =  , ( 12 ) where  is the azimuthal angle and r is the cylinder/circle radius. The normalized eigenfu nctions of 2 ,   are 41   1 , , , 1 , 0 , 1 , 2 in n en    =  = −  −  . ( 13 ) By plugging Eq. ( 13 ) into Eq . ( 12 ) we obtain 2 22 2 , 2 2 2 2 11 2 in n nn nn e r r r        − = = =  . ( 14 ) Thus, the eigenvalu es are given by (see Eq. (2) ) 2 2 n n r  = , which are used to buil d matrix Λ : ( ) 2 22 1 , , , , 1 , 0 , 1 , , , . r diag n n =   =   ΛΘ Θ ( 15 ) Note that the eig envalues are symm etric around the ground mode, 0 n = , an d that the ´spatial/angular ´ integrals of these eig enfunctions are zer o, except for t he ground mode: 8 ( ) ( ) ( ) ( ) 22 00 22 00 1 2 1 cos sin 2 2 , 0, 0, 0. in n d e d n d i n d n n               = =+  =  =       ( 16 ) For this geometry, the ele ments of matrix nm B (see Eq. (8); all the details are presented in Appendix A) are : ( ) ( ) 2 0 , 1 , 1 1 cos , 2 . 2 in im nm n m n m B r e e d r      − −+  =   =+  ( 17 ) This matrix can be r ewritten as 2 r = BK , ( 18 ) where K is a tri-diagonal matrix with 1 on the off-diagonals, , 1 , 1 1 n n n n KK −+ == , and 0 on the main diagonal, , 0 nn K = : 01 1 0 1 1 1 10      =       K . ( 19 ) 3.2 Optimizing fo r dimensio nality and speed: Re computing Λ and B In the previous section we derived the expressions for Λ and B for the basis of eigenfunctions. In practice, it is necessary to truncate the basis up to a maximum order nM = , resulting in matrices with dimensions o f ( ) ( ) 2 1 2 1 MM +  + , includin g the ground eigenmode. Since evaluating the matrix exponential terms in Eq. ( 11 ) is computationally demandin g when M is large, a more practical alternati ve is to build a ne w compact basis (i.e., n  ) with the same n umber of eigenfunctions nM = , but based on matrices with dimensions o f ( ) ( ) 11 MM +  + . This reduction is particularly advantageous from a computational point of v i ew. Modern alg o rith ms for evaluatin g the matrix exp onential, as implemen ted in MATLAB and SciPy , rely on a scaling-and-squaring strategy combined with Padé rational approximations and are domin ated by dens e matrix – matrix products, whose co mputational cost scales cubically with the matrix dimension 42,43 . Therefore, the p roposed reduction leads to an expected computational gain 9 proportional to ( ) ( ) ( ) 3 2 1 / 1 MM ++ . This corresponds to an ideal spe ed-up factor of eight in the asymptotic regi me of very larg e M . Th e following constructi o n ex p loits the symme try of the eigenfun ctions and the fact that the me asured signal is real. By adding the eigenfunctions with the sam e order n and opposite signs, we get the new functions ( ) ( ) 1 2 1 1 1 1 22 1 cos . n n n in in ee n        − − = + =+ = ( 20 ) Therefore, the new basis t o represen t the magnetizati on, as in Eq. (4), is defined as ( ) 1 , 0, 2 1 cos , 1. n n nn      =   =      ( 21 ) The normalization factor 12 added in Eq. ( 20 ) was chosen on purpose, to get an orthonor mal basis similar to Eq. (3) satis fying: ( ) ( ) 2 , 0 n m n m d        =  . ( 22 ) This is true in virtue of the following property: ( ) ( ) 2 0 2 , 0, cos cos , 0, 0, . nm n m d n m nm       ==   = =       ( 23 ) According to this operation, it is possible to transform the 21 M + vect or ( ) t c from the o rig inal basis to the compact basis c with 1 M + terms as : = c Tc , ( 24 ) where the ( ) ( ) 2 1 1 MM +  + transformation matrix is: 10 0 0 0 1 2 0 1 2 0 1 0 0 0 1 2 0 0 0 0 1 2        =        T . ( 25 ) By using this transfor mation, we rewri te Eq. ( 11 ) in terms of th e vector c as ( ) ( ) ( ) ( ) ( ) 0. D i G D D i G T t e e e      − − −  − − + = Λ B Λ Λ B c T Tc ( 26 ) Notably, in Appendix B, we prove (due to the structure of T and the symmetry of matric es Λ and B ) the following proper ty: 33 2 1 2 1 T TT T e e e e e e −− − − − − = A T A T A A T A T T A T TT . ( 27 ) Therefore, Eq. ( 26 ) can be rewritt en as: ( ) ( ) ( ) ( ) ( ) 0, D i G D i G D t e e e      − − − + −  − = Λ B Λ B Λ c c ( 28 ) where   ( ) ( ) 01 2 2 , , , , 1 , w h e r e 0 , 1 , , , , T diag d ia g n r  = =  = =   Λ T Λ T Θ Θ ( 29 ) and , , w he r e 2 02 2 0 1 1 1 10 T r = =      =       Β T ΒT K K ( 30 ) 11 The main benefit o f w orking with the c ompact basis, Eqs . ( 28 )-( 30 ), is that n ow the matrix e xponentials are e valuated on matrices with reduced dimensionality, i.e., ( ) ( ) 11 MM +  + instead of ( ) ( ) 2 1 2 1 MM +  + in the original basis, all owing for a m ore efficient i mplementation. 3.3 Diffusion con fined on cylind rical surfaces: Compu ti n g w and ( ) 0 c According to our pre vious results, th e dMRI signal E ⊥ is given by (see Eqs. ( 10 ), ( 28 )-( 30 )) ( ) ( ) ( ) 22 2 22 , , 0 D D G r D G r ii rr r E G e e e     −     − − − + −         ⊥   =      Θ K Θ K Θ wc , ( 31 ) where Θ and K are defined in Eqs. ( 29 ) and ( 30 ). However, to evaluate this expression we still need to determine w from the eigenfunctions (Eq. ( 21 )) , and ( ) 0 c from the initial c onditions. From Eq. ( 16 ) we see that the spatial int egrals in the diffusion ge ometry (on the circle ´s circumference) are zero 0 n w = , except for the element corresponding to 0 n = :   0 , 0 , w = w . This happens because when the ground eigenfunction ( ) 0  r is constant/spatially uniform, like in our problem, the spatial integrals of all the other eigenfu nctions are zer o due to the orthogonalit y property, s ee Eq. ( 22 ). On the other hand, vector ( ) 0 c can be determined fr om the initial conditions, that i s, fr o m the value of ( ) ,0 Mt = r . Assuming a spa tially uniform initial magnetization on the cylinder sur face implies ( ) 00 n c = for all 0 n  , i. e., ( )   0 1 , 0 , T = c . This follows from the expan sion in Eq. (4) for 0 t = , t ogether with the observation that the only spatiall y constant eigen function is the ground state ( ) 0  r . Thus, Eq. ( 31 ) becomes ( ) ( ) 22 2 22 0 0, 0 ,, D D G r D G r ii rr r E G w e e e     −     − − − + −         ⊥   =     Θ K Θ K Θ , ( 32 ) where we symbolically evaluated the vector multiplications. The subscript in the previous equation indicates that after computing the product of the three matrix exponentials we o nly take the entry at index ( ) 0, 0 , i.e., the first matrix element correspond ing to ground mode, the index-te rm 0 n = . Note that if we repea t the PGSE experiment for 0 G = , the above expressi on is simplifi ed to 12 ( ) ( ) ( ) 2 0 ,0 2 0 0 ,0 0 0 0, , , , , D r D r E G w e we w    + − ⊥ + −  =  =    = = Θ ( 33 ) where the m a trix-exponent ial b ecomes the ordinary exponential o f the r esulting diag onal matrix, and the first element of 0,0 0  = . This happens because in problems where ( ) 0  r is a constant, like in our case , 0 0  = is the only eigen value that fulfill s the equality given in Eq. (2), ( ) 2 , 0 0 0 0      = = r . Thus, Eq. ( 32 ) can be writte n in the more co mmon normaliz ed form: ( ) ( ) ( ) 2 2 2 22 0 ,0 ,, 0 D D qr D qr ii r r r Eq e e e Eq    − − + − − − ⊥ ⊥   =  =   Θ K Θ Θ K ( 34 ) which is dimensionl ess and the n ormalized sig nal lies in the unit inter v al [ 0,1 ]. N ote that we introdu ced the magnitude of the q-vector as qG  == q . This rela tionship provides an exact analytical expression for th e dMRI sig nal arisin g from diffusin g spin s in a cylindrical surface, when the gradient is applied perpendicularly to the cylinder´s axis . N otably, it is not based o n any simplifying assumptions about the diffusi on propagator, the duratio n of the gradient pulses, or the distribution of phases. The only approximation involved is numerical, result ing from truncating the number of eigenfunctions and eigenvalues employe d to practically evaluate th e expression (i.e., the dimensionality of matri ces Θ and K ). Neverthel ess, in practice, i t usually converge s with a fe w terms a s the eigenvalues gr ow as 2 n . 3.4 Exact reduction to a single matrix exp onential vi a conjugation symmet ry A direct evaluati on of Eq. ( 34 ) require s tw o matrix exp onentials p er sign al sampl e (e.g., g radient directi on or b -value). An additional reduction in computati onal cost follows from noting that the first and third matrix exponentials in Eq. ( 34 ) (corresponding t o the time intervals in which the diffusion g radients are applied in the P GS E experiment) are complex conju gates of each other. Therefore, in practice, it is sufficient to evaluate only one matrix exponential, while the second o n e can be obtained by complex conjugation. This strategy reduces the signal evaluation computational cost by ap proximately a factor of two. 3.5 Decoupled diffu sive motions an d dMRI signa l componen ts In the previ ous s ection we derived the analyti cal dMRI signal expr ession (Eq. ( 34 )) for the case when the diffusion-encoding grad ient is applied perpend icularly to the main cylinder´s axis. For an arbitrary gradient orientation, this e xpression is still valid but some m odifications mus t b e implemen ted. 13 For an infinitely long cylind er, the diffusion process o f spin particles parallel and perpendicular to the cylinder’s main axis is statistically independent. As a result, the dMRI signal can be expressed as the product of the signals ari sing from disp lacements par allel and perpendicu lar to t he cylinder’s axis 44 : ( ) ( ) ( ) , , , , , , E E E    ⊥⊥  =   q q q , ( 35 ) where q and ⊥ q are the projecti on of the q-space vector ˆ q = qg along the directions parallel and perpendicular to the cylind er’s axis, ⊥ =+ q q q ; thus their magnitudes are given b y ( ) c o s qq  == q and ( ) sin qq  ⊥ ⊥ == q , where  is the angle between the diffusion gradi ent o rien tat i o n unit vector ˆ g and the unit vector along the main cylinder’s ˆ v . A general detailed deriva tion of this type of decoupled signal model is pr ovided in 44 . Since the motion o f particles along the cylinder’s main axis is unrestricted, we assume 1D Gaussia n diffusion with a characteristic myelin water diffusivity D on the c ylinder’s surface. The normalized dMRI PGSE signal ( ) ,, E   q arising from th ese Gaus sian displacements is ( ) ( ) ( ) ( ) 2 2 3 cos ,, 0 qD E e E   −  −  = = q q , ( 36 ) where ( ) 2 3 bq  =  − is the b -value . On the other hand, in the expression for ( ) ,, E  ⊥⊥  q given by Eq. ( 34 ) w e must replace q by its projection q ⊥ along the plan e perpendicular to the cylin der axis. Therefore, the final signal expressi on in Eq. ( 35 ) for an arbitra ry g radient orientati on becomes ( ) ( ) ( ) ( ) ( ) ( ) ( ) 22 2 2 2 2 2 ˆˆ ˆˆ 11 ˆ ˆ 3 22 0,0 ˆ , , , 0 q r q r D D D ii qD r r r Eq e e e e Eq     −  −  − − + − − − −  −     =  =  g v g v Θ K Θ Θ K gv g , ( 37 ) where we substituted the following trigonometric relationships, ( ) ( ) 22 ˆ ˆ c o s  = gv and ( ) ( ) 2 ˆ ˆ sin 1  = −  gv , and made the signal to explicitly depend on the experimental (arguments) parameters of the PGSE sequ ence. 14 3.6 Strang splitt ing with p sub steps While in the previous sections we introduced theoretical developments and implementati o n aspects t o accelerate the ev alu ation o f the radial signal component, additional speed -u ps are desirable. This is relevant, for example, when a large number of signals m u st be evaluated repeatedly across diffusion gradient orientations, or across different radii when the signal model is used within a nonlinear fitting procedure (e.g., t o estimate the r adius from m easured data as an in verse proble m). In this section, we employ the s econd-order Strang splitting (or Strang – Marchu k operator spl itting) method 45,46 to obtain an approximate signal, which can significantly accelerate the signal evaluation while providing explicit con trol of the approximation accuracy. This approach is widely u sed in different fields, including the simulation of Hamiltonian systems, quantum mechanics (for approximating the time evolution governed b y the time-dependent Schrödi nger equation), as well as in fluid dynamics and numerical soluti ons of partial differential equations. The idea behind this method is to approximate the matrix exponential s in Eq. ( 37 ) as ( ) ( ) ( ) ( ) 2 2 22 2 ˆ ˆ 1 ˆ ˆ 1 2 2 2 2 2 , p qr D D qr D i i pr p pr r e e e Op e    − − − −  −    +    gv gv Θ Θ K ΘK ( 38 ) where the ap proximation error ( ) ( ) 2 Op  decreas es q uadratically with the number of subst eps p used to discretize the time interval  . Therefore, the err o r can be m ad e arbitraril y small by choosing a sufficiently large value of p . In our application, the computational gain arises from the fact that the m atric es Θ and K , which are fixed and do not depend on the radius or the sequence parameters (e.g., gradient orientation and grad ient strength), n ow app ear insi de short-time matrix exp onential te r m s. This allows us to diag o nali ze K (i.e., T = K Q ΩQ , where ( ) 0 , , M dia g   = Ω and T = Q Q I ) and write ( ) ( ) 22 ˆˆ ˆˆ 11 22 . qr qr ii T pp ee −  −   = g v g v K Ω QQ ( 39 ) Notably, this removes the need to evaluate a d ense matrix exponential; the exponential in the second term o f the pre vious equat ion is appli ed t o a diagonal m atrix an d can be c omputed by simply takin g the scalar exponen tial of each diagonal entry. Moreover, the diagonalization of K can be precompu ted before the signal ev aluation. By substituting Eq. ( 39 ) int o Eq. ( 38 ), we obtain ( ) ( ) ( ) 22 2 2 2 2 2 ˆˆ ˆˆ 11 ˆ ˆ 1 1 2 2 2 2 2 1 q r q r D D D qr D p ii i TT pr p pr p pr r i e e e e e e     −  −  − − − − −  −  =        g v g v gv Θ Θ Θ ΩΩ ΘK Q Q Q Q . ( 40 ) Combining Eq. ( 40 ) with Eq. ( 37 ), the PGSE signal can be efficiently evaluated. Notice that all m atric es involved in the resulting expression, except for the orthogonal matrix Q , are diagonal. In addition, the matrix products appearing in the signal expression do not need to be for med explicitly. Since the 15 measured signal corresponds to the first component of the resulting coefficient vector ( )   0 1 , 0 , T = c , it is sufficient to apply the successive matrix exponentials directly t o the initial v ec tor. As a result, all operations are carried out in a matrix – vect or form, which is co mputationall y cheaper than matrix – matrix prod ucts. This approximation was found to provide an accurate and substantially faster surrogate for repeated model e valuations f o r different gradien t ori entations. While the exact signal evaluation described in th e previous secti ons was used to generat e all re sults rep orted in the Results secti on, we includ e a dedica te d subsection to compare the approxi mate and exact signals as a function o f p , and to assess the corresponding c omputational speed-up , in order to evaluate the practical potential o f this appr oach for future data-driven app lications. 3.7 Spherical mean signal: Angular averaging and compu tational strategies In many d MRI ap plications, it is desirabl e t o work with observables that are invarian t to the orientation of the underlying microstructure. In this context, the spherical mean (powder-averaged ) dMRI signal — obtained b y averag ing th e signal over all gradient directi ons 47 – 49 — provides a rotationally invariant observable that depends only o n the intrinsic properties o f the confining geometry and the acquisition parameters. In the present work, we theref ore use the spherical mean sig nal to assess the beha vior o f the proposed model and to design efficient nu merical strategies for its e valuation. All angular average s reported in this w ork are comput ed using spherical Voronoi weights as sociated with the acq uisition g radient directions. Specifically, given a se t of unit gradient d irections   î g , the V oronoi- weighted spherical mean si gnal is computed as ( ) ii i E w E =  g , ( 41 ) where the weights i w correspond to the areas of the Voronoi cells o n the unit sphere associated with each gradient direction and are normalized such that 1 i i w =  . This weighted average mitigates the angular sampling b ias arising fro m no n-equid istant gradient direction s et s and provides a more accurate approximation of the continu ous spherical mean than uniform averag ing when the acquisition scheme is not perfectly unifor m. In addition to the scheme-averaged signal, we also compute the ideal spherical mean of the pr oposed model, defin ed as the cont inuous ang ular ave rage o v er the un it sph ere. Owin g t o the axial symmetr y of the cylindrical sur face mod el in Eq. ( 37 ), the spherical me an reduces t o a one-dimensional integral over ( )   co s 0 ,1  = : ( ) ( ) ( ) 22 22 2 2 2 11 1 3 22 0 0, 0 0 D q r q r D D ii qD r r r E e e e e d Eq      − −− − + − − − −  −   = =    Θ K Θ Θ K . ( 42 ) 16 This integral is evaluated by Gauss – Legendre quadratur e, which provid es a nu merically efficient and controllable appr oximation of the spherical mean. Unless otherwise stated, the exact analytical formulation derived in this work is used as numerical reference an d the spherical mean is computed using the Voron oi-weighted a v er age over the acquisition scheme . Accelerated im p lementations bas ed on th e Strang spl itting approach described in the previous subsection and on Gauss – Legendre angular quadrature for computing the spherical mean si gnal are evaluated with respe ct t o this referen ce in terms of accuracy and co mputational cost. 3.8 Monte Carlo s imulatio ns Monte Carl o Diffusion Si m u lations ( MCDS) were employed to validate the pr o p o se d analytical model. We used an MC simula to r developed by our group, available at https ://github.com/j onhrafe/Robust-Monte- Carlo-Simulations 50,51 . This tool has been validated ag ainst analytical models across multip le geometries, including imperm eable planes, cylin ders, spheres, an d more recently in cylind rical surfaces 35 . We compar ed the analytica l dMRI signal to those generated by the MC si m ulati o n s for identical cylindrical surfaces. Specifically, dMR I signals we re generated from 50 independent cylinders with radii uniformly spaced from 0.1 µ m to 5.0 µm in increments of 0.1 µ m following the protocol described in 35 . The diffusion process simulation was carried out u sing p N = 75 × 10 3 p articles uniformly distributed on each cylindrical surface and t N = 15 × 10 3 diffusion steps per particle. The axial myelin water diffusivi ty was set to D = 0.8 μm 2 /ms, as this valu e p roduced the larg est discrepanc ies in our previous model 35 . The dMRI da ta were gen erated using a PGSE sequence wi th trapezoidal diffusion g radients employin g a gradient strength of G = 500 mT/m and a sle w rate of SR = 500 T/m/s, which are c ompatible with the specifications of a Connectome 2. 0 3T scanner 52 an d small animal 7T and 9.4T Bio Spec MRI scanners. The protocol included 90° and 18 0° p ulse durations of 2 ms an d 4 ms , respectively . The data were g enerated for six b -values = [1, 2, 3, 4, 5, 6] ms / μm 2 , which wer e se lected u sing the shortest possible echo ti me ( TE ) for each case, while keeping the TE smaller than 20 ms an d main taining m a xim u m G and SR , following the implementation descri bed in 53,54 . For each b -value, dMRI signals were generated for 92 gradient orientations uniformly dist ributed on the unit sphere, along with the signal for b = 0. Table I shows the experimental para meters. The sub sequent analyses f ocused on the spherical mean sig nal normalized by the b = 0 signal. The sa me Voronoi weights ar e consistently applied to analytical m od el evaluations and Monte Carlo si mulations. Table I . Experimental p arameters for Monte Carlo simulatio n s using a P GSE sequence. b - value (ms/μm 2 )  (ms)  (ms) TE (ms) 1.0 7.72 2.88 13.43 2.0 8.72 3.89 15.44 3.0 9.45 4.62 16.90 4.0 10.03 5.20 18.06 5.0 10.53 5.70 19.06 6.0 10.97 6.14 19.94 17 4 Results In this secti on, we firs t illustrate charac te risti c feature s o f the diffusi on signal p redicted by the pr oposed cylindrical surface model and validate the exact analytical formulation against Monte Carlo simulations using scheme-averaged spherical mean signals comp uted with Voronoi weights. We then investig ate the numerical accuracy and c omputati onal efficiency of t he spec tral formulati on b y analyzing its convergence with respect to the trun cation order of th e spectral basis and by intr oducing two complementary acceleration strat egies based on opera tor splitting an d angular quad rature. 4.1 b-value – depend ent sphe rical mean signal profile Figure 1 illustrat es the theoretical spherical mean dMRI signal arisin g fro m diffu sion confined to a cylindrical surface, as predicted by the general analytical model introduced in Eq. ( 37 ). The signal was computed for a PGSE sequence with b -values ranging from 0 to 40 ms/µm², and three representati ve cylinder radii : 1.0 µm, 2. 0 µm, and 3. 0 µm. Additionall y, we plot the spherical mean dMRI signal from th e recently proposed spherical-mean Gau ssian phase approxi mation (GPA) model for this ge ometry 55 . Figure 1. Signal atten uation for d iffusion confined to cylindrical surfaces. The spherical mean dMRI signal was generated using the analytical m od el introduced in Eq . ( 37 ) f or a PGSE sequence with Δ = 1 1 ms and δ = 3 ms . The axial diffusivity was set to D = 0.8 µm²/ms. The signal attenuat ion is shown on a logarithmic scale for three c y li nder radii: 1 .0 µm (b lu e), 2 .0 µm (orange), and 3.0 µm (green). The analytical e xpression was eval u ated usin g M = 50 spectral m odes (continuous lin es) . Additionall y , we show the s pherical mean signal resulting f rom the Gaussian Phase Approximation (GPA) model 55 (dashed lines) . No te that th e large b -v alues used in this figure are intended to illustrate diffraction-like patter n s that are not necessar ily accessible with standard clinical proto cols. For small radii (e.g., r =1 .0 μm ), the GPA closely follows the exac t anal ytical sig nal over a wide range of b -values, with only minor deviations appearing at very high diffusion weightings ( b  20ms/μm 2 ). As the radius increases, however, noticeable differences between the two models emerge at lower b -values. For r =2 .0 μm and r =3 .0 μm , the GPA progressively devi ates from the exact soluti on for b  5 ms/μm 2 , 18 with th e discrepa ncy bec oming more pronounc ed f or larg er b -values. At very hig h diffu sion weig htings, the exact solution exhibits oscillatory features com monly ref erred to as diffraction p atterns, which ari se from the angular confinement of diffusion and th e resulting interference a mong discr et e Laplacian eigenmodes on the cylindrical surface. These features are not captured b y the GPA, highlighting the importance of the exact analytical for mulation when describing the signal in regimes where restricted diffusion effects b ecome dominant. 4. 2 Validat ion against Mo nte Carlo simulatio ns To validate the accuracy of the proposed analytical model, we compar ed its predictions against MC - generated signals . Analytical signals wer e computed u sing M = 50 modes. Figure 2 sho ws the spherical mean dMRI signal as a function of cylind er radiu s f o r an axial diff usivity of D = 0.8 µm²/ms, an d six b -values rang ing fr om 1.0 to 6.0 ms/µm². The analyti cal model (continuous lines) exhibits excellent agreement with the MC simulations (discrete points) across the entire ran ge of radii and b -values consider ed. This confirms the correctness of the spectral formulation an d the exac t treat ment of finite gradient pulse effects. As ex p ected, increasing the b -value leads to str o ng er s ignal attenuation and enh anced sensitivity to th e cylinder radius. For the acquisition protocol c onsidered here, the sign al d isplays li m ited sensitivity to radii smaller than approximately 0.5 µm and larger than approximately 4.0 µm, indicating a finite radius sensitivity windo w determi ned by the gradient streng th an d diffusion time s of the PGS E protocol. Figure 2. Spherical mean dMRI analytical versus MC signals. Spherical mean dMRI s ignals as a f unction of cy linder radius, co m puted usi n g the analytical model ( continu ous lin es) and Monte Car lo (MC) simulations (dots) are shown . Signals are evaluated for b -values of [1. 0, 2 .0, 3 .0, 4.0, 5.0, 6 .0] m s/µm², using a PGSE seq u ence with para m eter s listed in Table I and axial d iffu sivity D = 0.8 µm²/ms. Anal ytical signals a nd MC si m ulations w ere ge n erated f or 5 0 discrete radii spanning 0 – 5 µm. 19 4.3 Numerical an alyses To characterize the nume rical accurac y and compu tat ional performance of the p roposed analyti cal formulation, we investigat e three independent sources o f approximation error: (i) the trunca tio n order o f the spectral expansion, (ii) the nu mber of substeps used in the Strang splitting approximation, and (iii) the number of quadrature n o d es e mployed to e valuate the sph erical mean. In all cases, we quantify the associated accuracy-runtime trade-offs in order to assess the practical suitability of the proposed meth ods for repeated signal evaluati ons, as required i n parameter estimation and fitting applications. All experiments were perf ormed using the same PGSE parameters  = 10.97 ms and  = 6.14 ms, the same set of 92 diffusion-encodin g directions used in the Monte Carlo simulations, and the same three diffusion weightings b = [3, 6, 20] ms/μm². Unless otherwise stated, the spherical mean is computed using the Voronoi-weighted average over the gradient orientations. The analysis was carried out for r N =50 independent cylinder s with radii unifor mly distributed between 0.1 μm and 5.0 μ m . For all bench marks, exe cution was restricted t o a singl e CPU c ore and a sing le thread in order t o enable a fair comparison of the intrinsic computational cost, independent of hardware-specific parallelization or multi-threading effects. All computations were performed on a laptop equippe d with an Intel Core i7- 7700HQ CPU ( 8 cores, 2.80 GHz ) running Linux. 4.3.1 Convergen ce of the an alytical signal with respec t to th e truncation order We fir st investig ated the c onvergence of the anal ytical signal with r espect to th e truncation order M of the spectral expansi on. The Voronoi-weighted sph erical mean computed using M = 50 modes was ta ken as reference, denoted b y ref E . Voronoi-weighted spheri cal mean signals M E computed using lower truncation orders M = [3, 5, 10, 20] were compared ag ainst this referenc e. For each truncation ord er and each b -value, the total si gnal discrepancy over all r N =50 rad ii wa s quantified using the mean relative abs olute error ( MR AE ), defined as ( ) ( ) ( ) 1 1 100%. r N M i ref i M i r ref i E r E r MRAE N E r = − =  ( 43 ) In addition, we report the co mputation time r equired to g enerate the Vor onoi-weighted spherical m ean signal for all gradient dir ections per radius, together with the corresp onding accelerati on factor relati ve to the reference computati on using M =50 modes. For comparison, we also include the results obtained from the spherical- mean GPA m odel 55 .The resul ting accuracy and p erformanc e metrics ar e summarized in Table II. The results demonstrate very rapid convergence of the analytical signal with respect to the truncation order. In p articular, M =10 alr eady yields results that are numeri cally indis tinguishable from thos e obtained with sub stantially higher truncation o rd ers, while even M =5 pro vides acc urate appr oximations under the present acquisition conditions. This confirms that only a limited number o f angular m ode s 20 contribute significantl y to the signal, since the quad ratic growth of the Laplac e eigenvalues strongly suppresses higher- order contribu tions. For compari son, the spherical-mean GPA model provides the shortest comput ation time , yieldin g an acceleration factor of approximately 4.8×10 3 relative to the reference calculatio n. However, this gain in speed comes at the cost of substantially larger relative absolute errors, exceeding those o btain ed even with very low trunca tion orders (e.g ., M =3). Table II. Mea n relati ve absolute er ror ( MRAE ), computation time a nd acceleration factor as a function of t h e truncation order M of the spectral expansion. T he MRAE and acceleration factor are computed with respect to the Voronoi-weighted spherical mean signal obtained with M =50 m od es . Reported computation times correspond to th e average time required to ge n erate th e Voronoi-weig h ted spherical mean signal over the 92 dif f usion -encoding directions for a single cylinder radius and b -value. The average was computed b y dividin g the total runtime for all 50 radii and three diffusion weighting s by 50×3. R esults are repor te d for three diffusion weighti n gs: b = [3, 6, 20] ms/μm². Results f rom the spher ical mean Gaus s ian phase ap proximation (GPA) m o del are also included. Truncation order MRAE (100 %) Computation time ( milliseconds ) Acceleration factor (  ) b =3.0 ms/µm² b =6.0 ms/µm² b =20 ms/µm² M = 3 4.1 × 10⁻ 3 6.2 × 10⁻ 2 4.5 9. 8 ms 11.3 M = 5 2.7 × 10⁻ 7 2.4 × 10⁻ 5 5.1 × 10⁻ 2 10 . 2 ms 10.9 M = 10 3.8 × 10⁻ 11 1. 8 × 10⁻ 11 2 .5 × 10⁻ 10 12 . 5 ms 8.9 M = 20 3 .4 × 10⁻ 11 2.4 × 10⁻ 11 3.2 × 10⁻ 11 18 . 5 ms 6.0 M = 50 Reference Reference Reference 111 . 0 ms Reference GPA 3.05 14 .6 60 .0 0. 023 ms 4.8 × 10 3 4.3.2 Convergen ce of the ap proximated signal with r espect to the Strang splitting N ext , we assessed the convergence properties of the approxi mated analytical sig nal obtained using the Strang splitting app roach as a fun ction of the nu mber of substeps p . As reference, we used the Voronoi-weighted spherical mean signal generated with the exact analytical model using M=10 modes , which was selected bas ed on the convergence an alysis of the previous subsection. Approxima ted sign als comput ed using p =[5, 1 0, 20, 40, 80 ] substeps were compared against this reference. The anal ysis was perfo rmed for the same set of 50 cy linder rad ii and the sam e three b -values. F or eac h value of p , the overall d iscrep ancy with respect to the reference signal was quantified using the MRAE . In addition, we report the corresp onding computation times an d ac cele ra tio n f actors relative to the exact analytical evalua tion. The results ar e reported in Tabl e III. 21 Table III. Mean relative absolute er ror ( MRAE ), co m p utation ti m e and acceleration factor as a function of the num ber of Stra ng splittin g substeps p . The MRAE and acceleration factor are computed with respe ct to the exact analytical signal ge n erated w ith M =10 m od es. Repor ted computation times correspond to the average ti m e required to generate the Voronoi-weig h ted spherical m ean signal o v er the 92 diff u sion - encoding directions for a s ingle cylinder radius an d b -value. The average w as computed by dividi ng the total runtime for all 50 radii and three d if fusion weightings by 50×3. R esults are reported f or three diffusion weig h tings: b = [3, 6, 20 ] ms/μm². Splitting substeps MRAE (100 %) Computation time ( milliseconds ) Acceleration factor (  ) b =3.0 ms/µm² b =6.0 ms/µm² b =20 ms/µm² p = 5 0.55 1.01 2.0 0.6 ms 20 .8 p = 10 0.14 0.26 0.53 0.8 ms 15.6 p = 20 0.036 0.065 0.135 1.3 ms 9.6 p = 40 0.009 0.016 0.034 2.2 ms 5.7 p = 80 0.002 0.004 0.009 4.1 ms 3.0 The re sults in Table III sh o w that the Strang splitting approximation achieve s hig h accuracy with a relatively moderate numb er of substeps, while substantially reducing the computational cost. For example, for p =20 the MRAE is below 0.2% across all diffusion weightings, while the computation is approximately one order of magnitude faster than the exact evaluation. For p =40, the approximated signals are visually ind istinguish able from the exact sig nals (results not sho wn). As ex p ected, a s mall number of substeps leads to noticeable approximation errors, particularly at higher diffusion weigh tings, highligh ting the need f o r a minimal number o f substeps when opera ting in high- b regimes. 4.3.3 Convergen ce of the spherical mea n with respect to quadrature nod es Since the spherical mean constitutes a central observable in many applications, we investigated the convergence prope rties of the spherical mean sign al computed by Gauss – Legendre quad rature. All signals in this experime nt were generated using the exact analytical model with M =10 modes. The reference signal was taken as the scheme -a v erag ed signal computed using Voronoi weights for the 92 diffusion-encoding dir ections. This refer ence was comp ared with the spherical mean obtained by Gauss – Legendre quadratur e using q n =[5, 6, 7, 8, 20] nodes. For each value of q n and each b -value, the discrepancy with respect to the Vorono i-weighted spherical mean was quantified using the MRAE over all 50 radii employed in the previous sec ti o ns. The corresponding compu tation times and ac celeration fa cto rs are reported in Table IV. 22 Table IV. Mean relative absolute err or ( MRAE ), com putation time and acceleration factor as a fu nction of the number of quadrature nod es q n u sed in t he Gaus s – Legendre integration. All signals w er e ge n erated using t h e exact analytical model with M =10 modes. The MRAE and acce leration factor are computed with respect to the Voro n oi -weighted scheme-averaged signal ob tained u sing 92 diffusion - encoding directions. Reported com putation times co rrespond to the avera ge time required to generate the Voronoi -weighted spherical mean signal over the 9 2 diffusion -encoding directions for a single cylinder radius and b -value. The average was computed b y dividin g the total runtime for all 50 radii and three diffusion weighting s by 50×3. R esults are repor te d for three diffusion weightings: b = [3, 6, 20] ms/μm². Quadrature nodes MRAE (100 %) Computation time ( milliseconds ) Acceleration factor (  ) b =3.0 ms/µm² b =6.0 ms/µm² b =20 ms/µm² q n = 5 0.0262 0.0777 0.5867 0.7 ms 17 .9 q n = 6 0.026 2 0. 0776 0. 6282 0.9 ms 13.9 q n = 7 0.0262 0.0776 0. 6529 1.0 ms 12 .5 q n = 8 0.0262 0.0776 0. 6481 1.2 ms 10 .4 q n = 20 0.0 26 2 0.0776 0. 6485 2.8 ms 4.5 The r esults demonstrate stable and rapid con vergence of the Gauss – Legendre quadratu re with respect to the number of nodes. Very similar accuracy is observed for q n between 5 an d 20, indicating that a number of qu adrature nodes subst antially smaller than the number o f diffusion-encoding directions can be u sed in practice. This enables a significant reducti on in the computational cost of spherical mean evaluation, without requiring the explicit genera tio n of the signal for each experim ental gradient direc tion. After evaluatin g th e individ ual accele ration stra tegies in this and previous sections, we finall y assess their combined impact on the computation of the spherical mean signal. To this end, we performed a representative experiment in which the spherical mean an alytical signal wa s g enerated using a truncation order of M =10, evaluated through the Strang splitting approach with p =10 steps and q n =5 Gauss – Legendre quadrature nodes. Under these settings, the resulting MRAE values for the three considered diffusion w eightings ( i. e., b = [3, 6 , 20] ms/µm²) were [0.17, 0.34, 1.11] , respectively. The average computation time was 0 . 0 3 1 ms per cylinder radi us an d b -value. Th e implementation was further accelerated b y compiling the numeri cal routines usin g the just- in -time Numba c ompiler 56 . Comparing these results with those obt ained from the GPA model rep orted in Table I I for the same acquisition parameters shows that the c omputational cost of the accelerat ed analytical formulation is of the sa m e order of magn itude as that of th e GPA ( i. e., 0.031 ms vs 0.028 ms), while yielding substantially smaller approximati on errors. 23 5 Discussion The main result of this work is the exact analytical solution o f the Bloch – Torrey equation for diffusion confined to a cylindrical surface under finite PGSE gradients, yielding a closed-form expression for the corresponding dMRI signal valid for arbitrary gradient durations and separations. Notably, this solution does not rely on assu m p tions on the propagat or, spin phase distribu tion, or finite -p ulse effects. Related spectral approaches have been reported in the liter ature for restricted diffusion inside cylinders and spheres, and betw een parallel plan es with reflecting b o un daries 20,21,24,29,57 . The only appro ximation involved in the pra ctical evaluation o f the signal expression is nu merical, due t o the truncation of the sp ectral basis. This trun cation is well controlled, as the eig env alu es of the Laplac e operator on the cylindrical surface grow quad ratically with the mode index, lead ing to rapid convergence of the series. As a resul t, the analytical sig nal can be evaluated with high accurac y using a relatively small number of spectral mo d es (e.g., M =10; see Table II) . The proposed analytical signal expression was validated against Monte Carlo diffusion simulations over a wide range of cylinder radii and PGSE parameters (see Fig ure 2) . For the geometry under study, the diffusion process is co n strained along the (periodic) angular coordinate while remaining unbounded along the cylinder’s main axis. This leads to a natural deco up ling of the diffusive motion into perpe ndicular and parallel components, reflected in the factorized structure of the signal. We significantly reduced the dimensionality of the matrices involved in the signal evaluation by exploiting the symmetry of the cylindrical surface and introducing a reduced real spectral basis, thus achieving compu tational gain s without compromising accuracy. Beyond the exact analytic al formulation, an additional contribution of this work is the introducti on of efficient numerical strategi es that enable fast an d repeated evaluations of th e pr oposed signal model. In particular, we introduced a second-order Strang s plitting approxima tion for the matrix exponen tial operators ari sing in the PGSE sign al expression, which allows the tim e evolution to be decomp osed into a sequence of simple r op erators that can be efficiently appli ed thr ough a precomputed eigendecompositi on. This approach avoids repeat ed evaluations of dense matrix exponentials and enables the signal to be computed b y a sequenc e o f matrix – vector operations, yielding substantial computational savings while preser ving a controllab le approximation error (see Table III). In addition, we proposed an efficient evaluation of the spherical mean (powder-averaged) signal based on Gauss – Legendre quadrature, exploiting the axial symmetr y of the cylindrical surface model to reduce the angular integration to a one-dimensional quadrature problem (se e Table IV). T ogethe r, these developm ents provide a practical framework for accelera ting both the evaluati on of the directi onal signal and the computation of rotationally invariant observables, which are of central interest in many dMRI applications and in parameter estimation problems requiring rep eated model evaluati ons. The numerical analyses presented in this work further clarify the a cc urac y – efficien cy trade-offs associated with these approximation s. Although the results show that relatively small truncation orders o f the spectral expansion and a moderate number of Strang splitting substeps are sufficient to achieve high accuracy for the PGSE pro tocol considered here, these numerical parameters should not be adopted indiscriminatel y. In practic al ap plications, the requ ired truncation ord er depe nds on the acquisition parameters and on the range o f cylinder radii of interest, and a dedicated convergence analysis should therefore be carried out for each spe cific particular application. Likewise, the convergence o f the Stran g splitting approximation depends on the temporal structure and strength of the diffusion -en co din g 24 gradients, so that an application-specific assessment is needed to identify an appropriate numb er o f splitting substeps that balances accuracy and computational efficiency. Finally, the Voronoi-weighted powder average and the Gau ss – Legendre evaluation of the spherical mean correspond to two distinct numerical approximations of angular averaging. While the Voronoi-based approach provides a pr otocol- consistent estimate o f the powder-averaged signal for a finite and potentially non-uniform set of diffusion-encoding directions, the Gauss – Legend re quadratu re ap proximates the continuous, rotationally invariant spherical mean of the model. As a result, th ese tw o quantities do n ot necessarily coincide for a given acquisition scheme, and their differenc es reflect the fin ite and possibl y n on -uniform angular sampling of the pr otocol. In a previous study, we derived an exact analytical expression for diffusi o n constrained to cylindrical surfaces in th e narrow-puls e limit and pr o p osed an ap proximate extension to fin ite gradien t durations 35 . While that formulation captured the main features of surface-confined diffusion and showed good agreement with Mont e Carlo simulations in vari ous regimes, discrepancies emer ged for lar ge diffusi v itie s and high diffusion weightings, where finite-puls e a nd no n-Gau ssian effects become non-negligible. Recently, a new extension f or wid e pulses bas ed on the Gaussian phase approxim ation was in troduced in 55 , which y ield s accurate res ults for relatively small radii. The present work overcomes previous limitations by providing an exact treatment of finite gradient pulse durations, thereby completing the theoretical description present ed in previ ous works. From a broader theoretical perspective, a general spectral appr o ach was develope d to d escribe restricted diffusion in circular layers of arbitrary thickness under time -depend ent gradient fields 27 . Within such a framework, diffusion confined to a cylindrical surface can be formally interpr eted as a limiting case corresponding t o a vanishing radial layer thickness. T herefore, our work could be rederi ved following th e theoretical treatment pres ented in 27 . Ho wever, rather than relying o n this implicit limit, the present w ork follows a direct derivation tailored to the surface geometry. This strategy enables the explicit constructi on of the spectral geometry-specific m a trices governing the P GSE signal and yields a com pact analytical expression that is parti cularly suited for practical sign al evaluation. The present formulation relies on a number of idealizations that define its domain o f applicability. The cylinder is assumed to b e infinitely l o ng , thereby neglecting end effects. The surface is treated as p erfectly smooth and impermeable, and no exchange with other compartments is considered. In addition, relaxation effe cts and surface relaxi v it y are n o t explicit ly included in the model. While these assumpti ons are appropriate for isolatin g the effect of surfac e confinemen t on the diffusion sig nal, they should be kept in mind when considering applications to biological tissues or more complex geometrie s. The analytical framework developed her e can be extended in se veral directions. Possible extensions includ e the incorporation of distributions of (co nc entric) cylinder radii, alternative gradient waveforms, and additional physical effects s uch as relaxation or exchange. The presen ted model can b e extended to arbitrar y waveforms foll owing the app roach described in 20,21,30 . That is, by dividing the time interval of the waveform in to small subintervals and approximating the waveform by a piecewis e constant function in each subinterval. The optimal number of subinter vals depends on the smoothness o f the waveform. For example, a PGSE experimen t with trapezoidal gradients would require at least seven subintervals to represent the full wavefo rm, includi ng the rising and falling ramps of the tw o diffusion grad ients, the two plateau intervals at c onstant gradient ampli t u de, and the interval in which the gradient is zero. In practice, a fin er subdivision of the ramp in terv al s may be requ ired depending on the slew rate and on the desired accuracy of the waveform approxi mation. The resulting signal would involv e the multiplication of the matrix exponentials corresponding to each o f the 25 subintervals in which the waveform is assumed to be c onstant. As a result, although the prac tical implementation is s traightforward, the computation time grows p roportionall y with the nu mber o f subintervals. For this reason, we f o cused on rectangu lar PGSE gradi ents, as they involve three distinct subin tervals and allow an exact analytical treatment. In modern scanners with high slew rates, the ramp durations of trapezoidal grad ients are short c ompared to the pulse duration, so their impact o n the resultin g diffusion signal is rel atively small. As a result, trapezoidal gradi ents produce signals that are very well approximated by tho s e from rectangular gradients with the same no min al timing parameters ( ) ,   . Co n sistent with this, the signals generated by the proposed analytical mo d el were in excellent agreement with those obtained from Monte Carlo simulati ons using trap ezoidal gradien ts (see Table I an d Figure 2). While the cylindrical surface model is motivated by applications to myelin water diffusion 35,36,54,55,58 – 60 , this wor k f ocuses on the theor etical f ormulati on of the dMRI signal model, establishing a rigorous foundation for future experimental and data-driven studies. Closed- form analytical solutions of the Bloch – Torrey equation under fini te gradient pulses are availab le only for a limited number of geometri es. The present result extends t h is class to diffusi on confined to c ylindrical surfac es. Data availabili ty statement The dataset s and code supporting th e findings of th is study are publicly ava ilable at the following repository: https: //github.com/ejcan alesr/exact-PGSE- dMRI-signal-cylindrical -surface . Author contribution s EC -R: Conceptualizati o n, Investigation, Methodology, P r o jec t administration , Data curation, Forma l Analysis, Resources, S oftware, Supervisi on, Validation, Visuali zation, Writing – original draft, Writing – review and editing. CT: Inv estigation, Soft ware, Valid ation, Vi sualization, Writing – review an d editing . JM - G: In v estiga tio n , Resources, Writing – review and edi ting. DJ: Investigati o n, Resources, Writin g – review and editing. J-PT: Investigati on, Resources, Writing – review and editing. JR-P: In vestigation, Method ology, Data curation, Formal Analysis, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review and editin g. Acknowledgment s E.J.C.R. acknowledges finan cial support from the “Ramón y Cajal” Excellenc e Fellowship (Grant RYC2023- 042763-I ), funded by the Ministry of Science, Inn ovation and Un iversities (MICIU) and the Span ish Stat e Research Agency (AEI, 10.13039/50110001 103), and c o-financed by the Eur opean Social Fu nd Plus (ESF+ ). 26 Conflict of interest The authors declare tha t th e res earch was cond ucted without any commercial or fin ancial relationship s that could be constru ed as a potential c onflict of inter est. Generative AI sta tement The auth ors declare that C hatGPT- 5. 2 ( pa id version) wa s used to assist in identifying gram matical errors and typos in this manu script. All intelle ctual contributi o ns w ere made by the authors. Appendices Appendix A: Compu ting m atrix B The integral in Eq. ( 17 ) is so lved as: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 0 2 0 2 0 22 11 00 1 cos 2 1 cos 2 1 22 1 , 22 in im nm i n m i n m ii i n m i n m B r e e d r e d r e e e d r e d e d                  − − − − − − − +  =    =    +    =+       ( 44 ) where w e us ed Eul er´s for mula. W ithout l oss of gene rality, h ere w e assum ed that the diffusi o n gradien t is applied along the x-axis (assuming the cylinder ’s axis is parall el to the z-axis). Therefore, the angle between the gradient and t he position vector is  . Note this assumption does not affect the dMRI sig nal E ⊥ , sin ce it is constant fo r an y gradient orientation lying on the plan e p erpendicular to the cylinder’s axis . The above integrals can be evaluated by usin g the ide ntity: ( ) 2 , 0 2 i n m nm ed     − =  . ( 45 ) Thus, ( ) , 1 , 1 2 nm n m n m r B  −+ =+ . ( 46 ) Appendix B: Pro ving the relation ship in Eq. ( 27 ) The transformati on matrix in Eq. ( 25 ) satisfy the following properties: T = T T I is the identity and 27 1 / 2 0 0 0 1 / 2 0 1 / 2 0 1 / 2 0 0 0 1 0 0 0 1 / 2 0 1 / 2 0 1 / 2 0 0 0 1 / 2 T         =           TT . ( 47 ) Note that for a symmetri c matrix A , we get the following equality: ( ) T = AT T T AT . This can be demonstrated by multiplying both sid es of th e equati on by T T , from th e left, an d r eplacing T TT by the identity matrix. The relati o nshi p we want t o simpli fy is 3 21 T e e e − −− A AA TT . For this, let´s first write the matrix exp onential in Taylor series: ( ) 23 0 1 11 ! 2! 3 ! k k k e k  − = − = = − + − +  A A I A A A . ( 48 ) After substituting this relationshi p into the ter m 1 e − A we get ( ) ( ) ( ) 33 2 1 2 3 2 3 21 1 0 1 0 1 , ! 1 ! . T k T T k k k k TT k T e e e e e k ee k e e e  −− − − − =  − − = − −−  −  =    −  =   =   AA A A A A A A A T A T T T T A T T T T A T TT ( 49 ) In the previous deriva tion, we used the relati onship ( ) , k kT = A T T T AT ( 50 ) which can be obtained by doing the following iterative algebraic substitutions, ( ) ( ) ( ) ( ) 2 1 1 2 2 ... p k k T k T k T k p T − − − − − = = = = A A T A T T A T A A T T A T A T T AT A T T A T , until pk = . We can apply the sa me procedure to those matrix exponential s involving 2 A and 3 A , respectivel y, and Eq. ( 49 ) beco mes 33 2 1 2 1 3 21 3 21 . TT T TT T TT TT T e e e e e e e e e e e e −− − − − − − −− − −− = = = AA A A T A T T A T T A T T A T T A T T A T T A T T A T T T T T TT ( 51 ) 28 To o b tain the expression in Eq. ( 28 ), we sub stituted matrices ( ) 1 D i G  =+ A ΛB , ( ) 2 D  =  − A Λ , and ( ) 3 D i G  =− A ΛB in Eq. ( 51 ) with the matrices in the matri x exponenti al terms in Eq. ( 26 ), which are symmetric. References 1. Mitra PP, Sen PN. Effects of microgeo metry and surface relax ation on NMR pulse d-field- gradient experiments: Si mple pore geo metries. Phys Rev B . 199 2;45(1):143-156. doi:10.1103/PHYS REVB.45.14 3 2. Le Bihan D. The “wet mind”: water and functi onal neuroimag in g. Phys Med Biol . 2007;5 2(7):R57- 90. doi:S0031-91 55(07)80522-7 [pii]10. 1088/0031-9155 /52/7/R02 3. Brownstein KR, Tarr CE. I mportance of classical diffusion in NMR studies of water in biolo g ical cells. Phys Rev A . 1 979;19(6): 2446-2453. doi:10.11 03/PhysRevA.19.24 46 4. Basser PJ, Mattiello J, L eBihan D. MR diffusion tensor spectroscop y and imaging . Biophys J . 1994;66(1):25 9-267. doi:S000 6-3495(94)80775-1 [p ii]10.1016/S 0006-3495(94)807 75-1 5. Stejskal EO, Tanner JE. Spin Diffusion Measure m ent s: Spin Echoes in the Presence of a Time- Dependent Field Gradien t. J Chem Phys . 1965;42(1): 288. doi:10.106 3/1.1695690 6. Novikov DS, Kiselev VG, J espersen SN. On modeling. Magn Reson Med . 2018;79( 6):3172-3193. Accessed September 11, 201 8. http://doi.wiley.c om/10.1002 /mrm.271 01 7. Alexander DC, Dyrby TB, Ni lsson M, Zhang H. Imaging brain microstructure with diff usion MRI: practicality and appli cations. NMR Bio med . John Wiley & Sons, Ltd . 2019; 32(4):e3 841. doi:10.1002/nbm.3 841 8. Dyrby TB, Innocenti GM, B ech M, Lundell H. Validatio n S t rategies for the In terpretation of Microstructure I maging Usi ng Diffusion MRI . Vol 182. ; 2018:62-79. 9. Canales-Rodríguez EJ, Pizz olato M, Zhou FL, et al. Pore size esti mation in axon-mi micking microfibers with diffu sion-relaxation MRI. Magn Reson M ed . Published online 2024. doi:10.1002/MR M.29991 10. Barakovic M, Pizzola to M, Tax C MW, et al. Esti mating axon radius usin g diffusion-relaxation MRI: calibrating a surfac e-based relaxation model with histo logy. Front Neurosci . 2 023;17:1209521. doi:10.3389/FNIN S.2023.1209521 11. Veraart J, Nunes D, Rudrapatna U, et al. Noninvasi ve quantification of axon radii using diffusion MRI. Elife . 2020; 9. doi:10.7554/ eLife.49855 12. Jones DK, Alexander DC, Bowtell R, e t al. Microstructur al imaging of the human brain wi th a “super - scanner”: 1 0 key advantages o f ultra -strong gr adients for diffusion MRI. Neuroimage . Published online Ma y 22, 2018. doi :10.1016/j.neur oimage.20 18.05.047 13. Alexander DC, Hubbard PL, Hall MG, et al. Orientation ally invariant indices of axon diameter and density from diffusion MRI. Neuroimage . 2010;52(4):1 374-1389. doi:S105 3-8119(10)0077 5- 5 [pii]10.1016/j.neur oimage.20 10.05.043 29 14. Palombo M, Ianus A, Guerreri M, et al. SANDI: A c ompartment-based model f or non-invasive apparent soma and neurite imaging by diffusi o n MRI. Neuroimage . 2020;215. doi:10.1016/J.N EUROIMAG E.2020.116835 15. Jelescu IO, Palo mbo M, Bag nato F, Schilling KG. Chall enges for bi ophysical modeli ng of microstructure. J Neu r osci Methods . 2020 ;344:10886 1. doi:10.1016/J.JN EUMETH.2020.1 08861 16. Özarslan E, Shemesh N , Basser PJ, Ozarsl an E, Shemes h N, Basser PJ. A g eneral framewor k to quantify the effect of restricted diffu sion on the N MR signal with appli cations to double pu lsed field gradient NM R experiments. J Chem Phys . 2009; 130(10):10470 2. doi:10.10 63/1.3082078 17. Grebenkov DS. Pulsed-grad ient spin-echo monitoring of restricted diffusion in multilayered structures. J Magn Re son . 2010;2 05(2):181-195. doi :10.1016/j.jmr.2 010.04.017S10 90 - 7807(10)00119-9 [pii] 18. Nilsson M, Lasič S, Drobnjak I, et al. R esolution limit of cylinder diame ter estimati on by diffusion MRI: The impact of gradient wa veform and orientation disp ersion. NMR Biomed . 2017;3 0(7). doi:10.1002/nbm.3 711 19. Robertson B. Spin-echo d ecay of spins diffusin g in a bou nded region. Phys Re v . 1966;151(1) :273- 277. doi:10.110 3/PhysRev. 151.273 20. Caprihan A, Wang LZ , Fukushima E. A Multiple-Narrow-Pulse Approximati on for Restricted Diffusion in a Ti me-Varying Field Gradient. J Magn Res on Ser A . 1996;118( 1):94-102. doi:10.1006/JMRA. 1 9 96.0013 21. Callaghan PT. A simple matrix for malism for spin echo analysis of restricted diffu sion under generalized gradien t waveforms. J Magn Reson . 19 97;129(1):74-84. d oi:S1090-7807(97 )91233- 7 [pii]10.1006/jmre. 1997.1233 22. Callaghan PT, Codd S L. Generalised cal culation of NMR imaging edge effe cts arising fro m restricted diffusion in po r ous media. Magn Re son Imagin g . 1998;16(5-6):471-478 . doi:S0730- 725X(98)00070- 8 [pii] 23. Codd SL, Callaghan P T. Spin Echo Analysis of Restricted Diffusion under Generali zed Gradient Waveforms: Planar, Cylindrical, and Spherical Pores with Wall Relaxi v ity. J Magn Reson . 1999;137(2) :358-372. doi:10.100 6/jmre.1998.16 79S1090-7807(98)91679-2 [pii ] 24. Barzykin A V. The ory of Spin Echo in Restricted Geometries under a Step-wise Gr adient Pulse Sequence. J Magn Reson . 1999;13 9(2):342-353. 25. Axelrod S, Sen PN. N uclear magnetic res onance spin e choes for restri cted diffusion in an inhomogeneous field : Methods and as ymptotic regi mes. J Chem Phys . 2001; 114(15):687 8-6895. doi:10.1063/1. 1356010 26. Grebenkov DS. Multiexp o n ential attenuation of the C PMG spin echoes due to a g eometrical confinement. J Magn Reson . 2006 ;180(1):118-126. do i:10.1016/J.JMR. 2006.01.0 14 27. Grebenkov DS. Analy tical solution for restri cted diffusion in circu lar and sph erical layers und er inhomogeneous magnetic field s. J Chem Phys . 2008; 128(13):13470 2. doi:10.1063/ 1.2841367 28. Grebenkov DS. N MR survey of refl ected Brownian motion. Rev Mod Phys . 200 7;79(3):10 77-1137. doi:10.1103/REV MODPHYS.79.1077 /P1077_1_EQ 2.TIFF/MEDIU M 30 29. Grebenkov DS. Laplaci an eigenfun ctions in NMR. I. A n umerical tool. Concepts Magn Reson Part A . 2008;32A(4): 277-301. doi:10.10 02/cmr.a.20117 30. Grebenkov DS. Laplaci an eigenfun ctions in NMR. II. Th eoretical advances. Concep ts Magn Reson Part A . 2009;34A (5):264-296. doi:10.1 002/CMR.A.20 145 31. Lori NF, Conturo TE, Le Bihan D. Definiti o n of displacement probability and diffu sion time in q- space magnetic res onance measurements that us e fin ite-duration diffusion-enc oding gradients. J Magn Reson . 2003 ;165(2):18 5-195. doi:S10907807 03002805 [pii] 32. van Gelderen P, Despres D, van Zijl PC, et al. E valuation of restricted di ffusion in c ylinders. Phosphocreatine in rabb it leg muscle. J Magn Reson B . 1994;103(3):2 55-260. doi:10.1006/JMRB. 1994.1038 33. Balinov B, Jönsson B , Linse P, S öderman O. The N MR self-diffusion method appli ed to restric ted diffusion. Simulati on of echo attenuati on from molecul es in spheres and bet ween planes. J Magn Reson A . 1993;1 04:17-25. 34. Sö derman O, Jönsson B. R estricted diffusi on in cylindrical geo metry. J Magn Reso n . 1995;Series A,:94-97. 35. Canales-Rodríguez EJ, Ta x CMW, Fischi-Gomez E, Jones DK, Thiran JP, Rafael- Patiño J. A diffusion MRI model for rand om walks c onfined on cylindrical surfaces: towards n on-invasive quantification of myelin sheath radi us. Front Phys . 20 25;13:1516630. doi:10.3389/fph y.2025.1516630 36. Peled S, Cory DG, Raymond SA, Kirschner DA, Jolesz FA . Water Diffusion, T2, and Compartmentati on in Frog Sciatic Ner v e. Mag n Reson Med . 1 999;42(5):91 1. doi:10.1002/(sici) 1522-2594(19991 1)42:5<911 ::aid-mrm11>3.0.co; 2-j 37. Inouye H, Kirschne r DA. Evolution of m yelin ultrastruc ture and the major structural myelin proteins. Brain Res . 2 016;1641: 43-63. doi:10.1016 /J.BRAINRES.20 15.10.037 38. Inouye H, Kuo FH, Denninger AR, Wein hausen B, Burg ham m er M, Kirschner DA. Myelin structure in unfixed, single ner ve fibers: Scann ing X-ray microdi ffraction with a bea m size of 20 0 nm. J Struct Biol . 2017;2 00(3):22 9-243. doi:10.1016 /J.JSB.2017.0 7.001 39. Barzykin A. Exact soluti on of the T orrey-Bloch equatio n for a spin echo in re stricted geometries. Phys Rev B . 1998 ;58(21):14 171. doi:10.110 3/PhysRevB .58.14171 40. Torrey HC. Bloch Equ ations with Diffusion Terms. Phys Rev . 19 56;104(3):5 63-565. doi:10.1103/Ph ysRev.104.5 63 41. Griffiths DJ., Schroe ter DF. Introduction to Quantum Mechanics . Third edit. Cam bridge Universi ty Press; 2018. https://www.librar ysearch.manches ter.ac.uk/discov ery/fulldisplay ?vid=44MAN_I NST:MU_NUI&t ab=local&docid=alma 992979 324700501631&lang =en&context=L 42. Higham NJ. The Scaling and Squaring Method for the Matrix Exponential Revisited . https://doi.org /101137/0406 1101X . 2006;26(4):117 9-1193. doi:10.11 37/04061101X 43. Al -Mohy AH, Higham NJ . A New Scaling and Squaring Algorithm for the Matrix Exp onential. https://doi.org /101137/0907 4721X . 2009;31(3):970-9 89. doi:10.1137 /09074721X 31 44. Assaf Y, Freidlin RZ, R ohde GK, Basser PJ. New modeli ng and experimental framework t o ch aracterize hin dered and restrict ed water diffusi on in brain white matter. Magn Reson Med . 2004;52(5):96 5-978. doi:10.1002 /mrm.20274 45. Strang G. On the C onstruction and C omparison of Diffe rence Schemes. https://doi.org /101137/0705 041 . 2006;5(3):506-517. doi:10.1137/0 705041 46. Marchuk GI. Some application of spli tting-up methods to the solution of mathematical ph ysics problems. Apl Ma t . 13(2):103-132. 47. Edén M. Computer simulations in s olid-state NMR. III. Powder averaging . Concepts Magn R eson Part A . 2003;18A (1):24-55. doi:10.1 002/CMR.A.1006 5 48. Lasič S, Szczepankiewi cz F, Erikss on S, Nilsson M, T opgaard D. Microaniso tropy im aging: quantification of microscopic di ffusion anisotropy and o rienta tio nal o rder parameter by diffusion MRI with magic-ang le spinning of the q-vector. Front Phys . 20 14;2(February). doi:10.3389/fph y.2014.00011 49. Kaden E, Kelm N D, Carson RP, D oes MD, Alexander DC. Multi-compart ment microscopic diffu sion imaging. Neuroimage . 2016;139: 346-359. doi:10.10 16/j.neuroi mage.2016.06.0 02 50. Rafael-Patino J, R omascano D, Ra mirez-Manzanares A, Can ales-Rodríguez EJ, Gi rard G, Thiran JP. Robust Monte-Carl o Simulations in Diffusion-MRI: Eff ect of the Substrat e Complexity and Parameter Ch oice on the Reproducib ility of Result s. Front Neuroinform . 20 20;14:457195. doi:10.3389/FN INF.2020.00008 /BIBTEX 51. Villarreal-Haro JL, Gardier R, Canales-Rodríguez EJ, et al. CACTUS: a co mputational framew ork for generating realis tic white matter microstructure subs trates. Front Neuroin form . 2023;17:1208 073. doi:10.3389/FNINF. 2023.1208073 /BIBTEX 52. Huang SY, Witzel T, Keil B, et al. C onnectome 2. 0: Developing the next-generation ultra- high gradient strength human MRI scann er for bridging studies of the micro-, meso- and macro- connectome. Neuroimag e . 2021;243. d oi:10.1016/J.NEUROIMA GE.2021.118 530 53. Mueller L, Tax CMW , Jones DK. Unprecedented echo times for diffusion MRI usin g connectom gradients, spiral r eadouts and field monitoring. MAGN ETOM Flash . Published online 2 019. Accessed December 14, 2023. https://cdn0.scrvt.co m/39b415fb 07de4d9656c7b5 16d8e2d907/ 18000000062 77360/3a1f90ead0 87/siemens-healthin eers-magnetom-flash - 74 -ismrm-unprecedented-e cho- times_180000 0006277360.p df 54. Tax CMW, Rudrapatna US , Mueller L, J ones DK. Chara cterizing diffusion of myelin water in th e living human brain usin g ultra-strong gradi ents and spiral readout. In: ISMRM 27th Ann ual Meeting & Exhib ition, 11-16 May 201 9, Montréal, QC, Canada. Abstract # 1115 . 2019. Accessed July 19, 2024. ht tps://archive.is mrm.org/201 9/1115.html 55. Lee HH, Chan KS, No vikov DS, Fieremans E, Huang SY. Axon Diameter Map ping from Myelin Water Diffusion M RI. IEEE Trans Med I maging . Publis hed online 2026. doi:10.1109/T MI.2026.36643 28 56. Lam SK, Pitrou A, S eibert S. Numba: A LLVM-based Pyth on JIT Compiler. Proc LLVM- HP C 2015 2nd Work LLVM Compil Infrastruct HPC - Held con junction with SC 2 015 Int Conf High Perform Comput Networking, Storage An al . 2015;201 5-January. doi:10. 1145/28331 57.2833162 32 57. Sukstanskii AL, Yabl onskiy DA. Effe cts of restricted dif fusion on MR signal formation. J Magn Reson . 2002;157( 1):92-105. doi:S10 907807029258 26 [pii] 58. Avram A V., Guid on A, Song AW. Myelin water weight ed diffusion tensor imaging. Neuroimage . 2010;53(1):13 2-138. doi:10.1016 /j.neuroimage.20 10.06.019 59. Brusini L, Menegaz G , Nilsson M. Monte Carlo Simula tions of Wa ter Exchange Thr o u gh Myelin Wraps: Implicati ons for Diffu sion MRI. IEEE Trans Med Imaging . 2019;3 8(6):1438-1445. doi:10.1109/T MI.2019.28943 98 60. Andrews TJ, Osborn e MT, Does MD. Diffu sion of myelin water. Mag n Reson Med . 2006;56(2):38 1-385. doi:10.1002 /MRM.20945