= 1mm

One-loop corrections to the spectral action

Teun D.H. van Nuland

Walter D. van Suijlekom

Institute for Mathematics, Astrophysics and Particle Physics

Radboud University Nijmegen, Heyendaalseweg 135,

6525 AJ Nijmegen, The Netherlands

We analyze the perturbative quantization of the spectral action in noncommutative geometry and establish its one-loop renormalizability in a generalized sense, while staying within the spectral framework of noncommutative geometry. Our result is based on the perturbative expansion of the spectral action in terms of higher Yang–Mills and Chern–Simons forms. In the spirit of random noncommutative geometries, we consider the path integral over matrix fluctuations around a fixed noncommutative gauge background and show that the corresponding one-loop counterterms are of the same form so that they can be safely subtracted from the spectral action. A crucial role will be played by the appropriate Ward identities, allowing for a fully spectral formulation of the quantum theory at one loop.

§ INTRODUCTION

Noncommutative geometry [15] offers a spectral viewpoint to geometry that allows to simultaneously capture field theories and gravity in a single framework. In fact, it allows for a unified geometrical derivation of the Standard Model of particle physics minimally coupled to gravity [10, 39], including the Higgs mechanism and the see-saw mechanism to yield masses for the right-handed neutrinos. This extends beyond the Standard Model to yield Pati–Salam grand unification [12, 13], which is currently one of the few candidate BSM-theories that is still found to be compatible with experiment. Variations on particle theories obtained in the same framework are considered in [8, 22, 21, 23, 7, 4, 20, 19, 35, 6], while the more foundational aspects on quanta of geometry were considered in [11].

The key ingredient in this description of field theories arising from noncommutative spaces is the spectral action principle [9]. It yields Lagrangians that are based solely on the spectrum of a given Dirac operator on a noncommutative spacetime. In the applications to particle physics phenomenology one then adopts the usual renormalization group methods to arrive at couplings and mass parameters at lower energy. Even though the appearance of such experimentally testable results from a geometrical framework valid at high-energies is very intriguing, we must confess that this step is a weak point of the noncommutative approach to particle physics. Indeed, it means that in the passage to the quantum theory one looses the elegant spectral and unifying picture that one started with and which one admired so much.

In this paper, we take a crucial step in the quantization program and analyze the form of loop corrections to the spectral action. Working in a very general context, in fact beyond [18, 14], we find that the resulting quantum fluctuations can be entirely formulated within the same unifying spectral framework and is thus a major improvement with respect to the usual RG-approach to the spectral action. The approach we take to the perturbative quantization of the spectral action is that  of random noncommutative geometries [25, 2, 32] (see also [3, 24] for computer simulations). More specifically, we adopt the background field method for which the path integral will be defined over all matrix fluctuations around a fixed noncommutative gauge background.

The key mathematical input is given by our paper [40], which gives a perturbative expansion of the spectral action in terms of noncommutative integrals over higher Yang–Mills and Chern–Simons forms.
We will here show that the one-loop corrections to the spectral action are of exactly the same form, and can thus safely be subtracted as counterterms from the spectral action. This establishes one-loop renormalizablity in the generalized sense of , where one allows for infinitely many counterterms.

§ DIAGRAMMATIC EXPANSION OF THE SPECTRAL ACTION

The spectral action [9] is defined on the eigenvalue spectrum $\{\lambda_k\}_k$ of a Dirac operator $\dirac$ by

\tr f(\dirac) = \sum_{k} f(\lambda_k)

for some suitable even function $f$. We want to analyze the spectral action for perturbations $\dirac \to \dirac+V $ by bosonic gauge fields of the form $V =a_j \diraccom{b_j}$ (summation over $j$ understood), where $a_j,b_j$ are coordinate functions on a noncommutative space. Even though our analysis is valid in the general setting of noncommutative geometry [15, 16] the most interesting cases that occur in physics are:

* Hermitian matrix models where both $\dirac$ and $V$ are hermitian matrices.

* Almost-commutative geometries $M \times F$, where $M$ is the spacetime manifold with Dirac operator $\slashed{\partial}$ and $F$ is a discrete noncommutative space describing the internal degrees of freedom, also equipped with a `finite' Dirac operator $D_F$. The gauge fields $V$ describe both Yang–Mills gauge fields $\mathbb{A}$ and scalar (Higgs)  fields $\Phi$ in the sense that

V= \underbrace{ a_j \slashed \partial (b_j)}_{\slashed{\mathbb{A}}} + \gamma_5 \underbrace{ a_j [D_F,b_j]}_{\Phi},

More details, also on the applications to particle physics, can be found in [10, 12, 22, 39, 14].

Our starting point is the following expansion of the spectral action [37, 36, 40]:

\begin{equation}
\label{eq:exp}
\tr\left( f(\dirac+V) - f(\dirac) \right) = \sum_{n=1}^\infty \frac 1 n \langle \underbrace{V,\ldots,V}_{n}\rangle.
\end{equation}

The brackets stand for the following contour integrals:

\langle V_1, \ldots, V_n \rangle \\=  \tr  \oint \frac{dz}{2\pi i} f'(z)  V_1 (z-\dirac)^{-1}   \cdots V_n (z-\dirac)^{-1}

where $V_1,\ldots, V_n$ are gauge fields as above; this can be represented nicely as a Feynman diagram:

\begin{equation}
\langle V_1, \ldots, V_n \rangle =
\scalebox{.7}{ \parbox{6cm}{
\hspace{1cm}
\begin{fmfgraph*}(30,30)
\fmfset{arrow_len}{3mm}
\fmfsurroundn{e}{6}
\fmflabel{$V_1$}{e4}
\fmflabel{$V_2$}{e3}
\fmflabel{$V_3$}{e2}
\fmflabel{$V_4$}{e1}
\fmflabel{$V_n$}{e5}
\fmf{photon}{e1,i1}
\fmf{photon}{e2,i2}
\fmf{photon}{e3,i3}
\fmf{photon}{e4,i4}
\fmf{phantom}{e6,i6}
\fmf{photon}{e5,i5}
\fmfrcyclen{fermion,left=0.25,tension=2}{i}{6}
\fmf{phantom,left=.5,tension=0.3}{e1,e5}
\fmfipath{p}
\fmfiset{p}{vpath(__e1,__e5)}
\fmfi{dots}{subpath (length(p)/4,3*length(p)/4) of p}
%        \fmf{dots,left=.5,tension=0.3}{e1,e5}
\end{fmfgraph*}}}
\label{eq:bracket}
\end{equation}

The loop diagram nicely reflects the cyclicity of the bracket: $\br{V_1,\ldots,V_n}=\br{V_n,V_1,\ldots,V_{n-1}}$. The second crucial property is that

\begin{equation*}
\br{a V_1,\ldots,,V_n}-\br{V_1,\ldots,V_n a}
=\br{\diraccom{a}, V_1,\ldots,V_n}
\end{equation*}

for any (noncommutative) coordinate function $a$.
In fact, this identity boils down to the following Ward identity,

(z-\dirac)^{-1} a - a (z-\dirac)^{-1} = (z-\dirac)^{-1} \diraccom {a} (z-\dirac)^{-1},

and may be represented diagrammatically:

\begin{equation}
\label{eq:ward}
\parbox{1.5cm}{    \begin{fmfgraph*}(10,10)
\fmfset{arrow_len}{3mm}
\fmfleft{lb,l}
\fmfright{rb,r}
\fmf{fermion}{l,r}
\fmf{dashes}{r,rb}
\fmfv{label=$a$,label.angle=-90}{rb}
\end{fmfgraph*}}
- \quad
\parbox{1.5cm}{       \begin{fmfgraph*}(10,10)
\fmfset{arrow_len}{3mm}
\fmfright{rb,r}
\fmfleft{lb,l}
\fmf{fermion}{l,r}
\fmf{dashes}{lb,l}
\fmfv{label=$a$,label.angle=-90}{lb}
\end{fmfgraph*}}=
\parbox{1.5cm}{
\begin{fmfgraph*}(15,10)
\fmfset{arrow_len}{3mm}
\fmfleft{lb,l}
\fmfright{rb,r}
\fmfbottom{cb}
\fmf{fermion}{l,c,r}
\fmffreeze
\fmf{photon}{cb,c}
\fmflabel{$\diraccom{a} $}{cb}
\end{fmfgraph*}}
\end{equation}

§.§ The brackets as noncommutative integrals

We want to express the amplitudes corresponding to the above loop diagrams in terms of suitable noncommutative integrals [15] (cf.  <cit.> or <cit.>. They are defined by

\begin{align}
\label{eq:bracket-cochain}
\int_{\phi_n} a^0 da^1 \cdots da^n %&= \br{a_0 D(a_1) , D(a_2),\ldots, D(a_n) }\\
&:= \scalebox{.7}{ \parbox{6cm}{
\hspace{1.5cm}
\begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{3mm}
\fmfsurroundn{e}{6}
\fmflabel{$a^0\diraccom{ a^1}$}{e4}
\fmflabel{$\diraccom{ a^2}$}{e3}
\fmflabel{$\diraccom{ a^3}$}{e2}
\fmflabel{$\diraccom{ a^4}$}{e1}
\fmflabel{$\diraccom{ a^n}$}{e5}
\fmf{photon}{e1,i1}
\fmf{photon}{e2,i2}
\fmf{photon}{e3,i3}
\fmf{photon}{e4,i4}
\fmf{phantom}{e6,i6}
\fmf{photon}{e5,i5}
\fmfrcyclen{fermion,left=0.25,tension=2}{i}{6}
\fmf{phantom,left=.5,tension=0.3}{e1,e5}
\fmfipath{p}
\fmfiset{p}{vpath(__e1,__e5)}
\fmfi{dots}{subpath (length(p)/4,3*length(p)/4) of p}
\end{fmfgraph*}}}.
%  \parbox{5cm}{\includegraphics[scale=.5]{ncintegral.png}}
\end{align}

For example, for one external edge we find

\begin{align}
\br{V} =   \br{a_j \diraccom{b_j}} &=%
\scalebox{.7}{  \parbox{2  cm}{
\begin{fmfgraph*}(20,20)
\fmfset{arrow_len}{3mm}
\fmfleft{l}
\fmfright{r}
\fmf{photon,tension=.5}{l,c}
\fmfv{label=$a_j\diraccom{b_j}$,label.angle=-90}{l}
\fmf{phantom,tension=.1}{r,c}
\fmf{fermion,left,tension=.5}{r,c,r}
\end{fmfgraph*}}}
= \int_{\phi_1} A,
\end{align}

where we have defined $A= a_j d b_j$ as the universal gauge form underlying the physical gauge field $V = a_j \diraccom{b_j}$. Note that the vanishing of this tadpole diagram corresponds to the vanishing of the first derivation of the spectral action under perturbations $D \mapsto D+V$. For natural choices of $D$ one may thus expect this term to vanish and, in fact, [18] works under this `vanishing tadpole' assumption.

For two external edges, we apply the Ward identity (<ref>) and derive

\begin{align*}
\br{V,V}&=
%\br{a_j \diraccom{b_j},a_{j'} \diraccom{b_{j'}}} =
\scalebox{.7}{\parbox{6cm}{
\hspace{1.3cm}
\begin{fmfgraph*}(15,15)
\fmfleft{l}
\fmfright{r1,r2}
\fmf{photon}{l,c1}
\fmfset{arrow_len}{3mm}
\fmf{fermion,left=1,tension=.5}{c1,c2,c1}
\fmf{dashes}{r2,c2}
\fmf{photon}{r1,c2}
\fmfv{label=$a_j\diraccom{ b_j}$,label.angle=-120}{l}
\fmflabel{$\diraccom{ b_{j'}}$}{r1}
\fmflabel{$a_{j'}$}{r2}
\end{fmfgraph*}} }\\[8mm]
&=\qquad \scalebox{.7}{\parbox{2.7cm}{
\begin{fmfgraph*}(15,15)
\fmfright{r}
\fmfleft{l1,l2}
\fmf{photon}{r,c1}
\fmfset{arrow_len}{3mm}
\fmf{fermion,left=1,tension=.5}{c1,c2,c1}
\fmf{dashes}{l2,c2}
\fmf{photon}{l1,c2}
\fmflabel{$\diraccom{ b_{j'}}$}{r}
\fmflabel{$a_j\diraccom{ b_j}$}{l1}
\fmflabel{$a_{j'}$}{l2}
\end{fmfgraph*}}}
+  \qquad \scalebox{.7}{\parbox{3cm}{
%    \hspace{1.5cm}
\begin{fmfgraph*}(15,15)
\fmfsurroundn{e}{3}
\fmflabel{$a_j\diraccom{ b_j}$}{e3}
\fmflabel{$\diraccom{ a_{j'}}$}{e2}
\fmflabel{$\diraccom{ b_{j'}}$}{e1}
\fmf{photon}{e1,i1}
\fmf{photon}{e2,i2}
\fmf{photon}{e3,i3}
\fmfset{arrow_len}{3mm}
\fmfrcyclen{fermion,left=0.5,tension=.5}{i}{3}
\end{fmfgraph*}}}
\\[6mm]
&= \int_{\phi_2} A^2 %a db  a' db'
+ \int_{\phi_3} A d A.% a db  d a' db'
\end{align*}

Similarly, by applying the Ward identity several times one finds that [40]

\begin{align*}
\br{V,V,V}&=\int_{\phi_3}A^3+\int_{\phi_4}AdAA +  \cdots %+\int_{\phi_5}AdAdA
\br{V,V,V,V}&=\int_{\phi_4}A^4+ \cdots . %\int_{\phi_5}(A^3dA+AdAA^2)+\int_{\phi_6}AdAdAA+\int_{\phi_7}AdAdAdA
\end{align*}

We now introduce a noncommutative integral $\int_\psi$ that differs from $\int_\phi$ by a total derivative:

\begin{equation}
\label{eq:psi}
\int_{\psi_{2k-1}}  \omega =
\int_{\phi_{2k-1}}  \omega - \frac 1 2 \int_{\phi_{2k}}  d \omega
\end{equation}

and rewrite the above brackets in terms of $\psi_1$ and $\psi_3$, as well as the remaining $\phi_2$ and $\phi_4$.
For the first two terms, we readily find

\int_{\phi_1}A + \frac 12 \int_{\phi_2}A^2 =
\int_{\psi_1}A + \frac 12\int_{\phi_2} (dA +A^2).

while after a slightly more involved derivation we also find for the next few terms that

\begin{align*}
&  \frac 12 \int_{\phi_3}AdA+ \frac 1 3 \int_{\phi_3}A^3+ \frac 1 3 \int_{\phi_4}AdAA + \frac 1 4 \int_{\phi_4}A^4
\\
% &=\frac 12 \int_{\psi_3} \left ( AdA+ \frac 2 3 A^3 \right) + \frac 14 \int_{\phi_4} dA dA + \frac 16 \int_{\phi_4} d(A^3) \\
% &\qquad   + \frac 1 3 \int_{\phi_4}AdAA + \frac 1 4 \int_{\phi_4}A^4\\
&\quad= \frac 12 \int_{\psi_3} \left ( AdA+ \frac 2 3 A^3 \right) + \frac 14 \int_{\phi_4} (dA+A^2)^2+ \cdots .
\end{align*}

The important message from the above derivation is that the expansion of the spectral action yields Yang–Mills and Chern–Simons terms. In fact, if we write $F = dA+A^2$ for the curvature and

\begin{gather*}
\cs_1(A) = A; \qquad  \cs_3(A) = \frac 12 \left( A dA + \frac 2 3 A^3 \right),
%\cs_5(A) = \frac 13 \left( A (dA)^2 + \frac 3 4 A dA A^ 2 + \frac 3 4 A^3 dA + \frac 3 5 A^5 \right),
\end{gather*}

then it turns out that the expansion has the following form of a Yang–Mills–Chern–Simons theory:

\begin{align*}
S_D[V] &= \int_{\psi_1} \cs_1(A) + \frac 12 \int_{\phi_2} F \\
&\qquad + \int_{\psi_3} \cs_3(A) +\frac 1 4 \int_{\phi_4} F^2% \\
%&\qquad + \int_{\psi_5} \cs_5(A) +\frac 1 6 \int_{\phi_6} F^3
+ \cdots
\end{align*}

Quite surprisingly, the systematics behind this derivation persists at all orders [40], while being based solely on the cyclicity of the loop diagram and the Ward identity (<ref>). It yields the following expansion for the spectral action

\begin{equation}
\label{eq:YmCs}
S_D[V] = \sum_{k=1}^\infty \left( \int_{\psi_{2k-1}}  \cs_{2k-1} (A) +\frac 1 {2k} \int_{\phi_{2k}}  F^{k} \right).
\end{equation}

The higher-order Chern–Simons forms are defined as in <cit.> by

\begin{equation}\label{eq:cs}
\cs_{2k-1}(A) := \int_0^1 A (t dA + t^ 2 A^2)^{k-1} dt.
\end{equation}

Again based solely on cyclicity of the loop diagram and the Ward identity, one can show that the integrals over $\phi_{2k}$ and $\psi_{2k-1}$ define even and odd cyclic cocycles, respectively; we refer to [40] for more details.

§ LOOP CORRECTIONS TO THE SPECTRAL ACTION

In order to analyze the quantum theory corresponding to the above classical action functional $S_D[V]$ we adopt the background field method. That is to say, we take the background fields to be gauge fields of the form $V =a_j \diraccom{b_j}$. However, the path integral is defined over the ensemble of all finite-size hermitian complex-valued matrices. This is in the spirit of random noncommutative geometries in the sense of [25, 2, 32] (see [3, 24] for computer simulations). As in these works, we consider the dimension, say $N$, of these matrices as a regularizing cutoff of our model, which should eventually be sent to $\infty$, while allowing us to realize our quantum theory as a hermitian matrix model.

In fact, for such finite-size matrices $Q =(Q_{kl})$, the brackets can be conveniently expressed in terms of divided differences of $f'$ [37]:

\begin{align*}
\frac 12   \langle Q, Q\rangle &= \frac 12\sum_{k,l} Q_{kl} Q_{lk} f'[\lambda_k,\lambda_l]\\
\frac 13   \langle Q, Q,Q \rangle &= \frac 13 \sum_{k,l,m} Q_{kl} Q_{lm} Q_{mk} f'[\lambda_k,\lambda_l, \lambda_m]
\end{align*}

et cetera, where $\lambda_k$ are the eigenvalues of $D$. Recall that the first two divided differences are defined by $f'[x,y] = (f'(x)-f'(y))/(x-y)$ and $f'[x,y,z]= (f'[x,y]-f'[y,z])/(x-z)$.

We now make the assumption that the first divided difference of $f'$ is strictly positive on the $N$ relevant eigenvalues of $D$ (see Figure <ref>). We may then perform the Gaussian integration as in, without the need for introducing a gauge-fixing and ghost sector, to get for the propagator:

\wick{\c Q_{kl} \c Q_{mn} }=\frac{ \int Q_{kl}  Q_{mn} e^{-\frac 12 \br{Q,Q}} dQ} {\int e^{-\frac 12 \br{Q,Q} } dQ} =
\delta_{kn}\delta_{lm} G_{kl}

in terms of $G_{kl} := \frac 1 {f'[\lambda_k, \lambda_l]}$. Notice that the inverse propagator is bounded, which is in stark contrast to the usual unbounded nature of inverse propagators in ordinary local quantum field theory. We see this as another manifestation of the regularizing properties of the spectral action, in line with [38, 29, 33, 1].

It is an interesting problem to analyze the form of the propagator for more general $f$, including a possible gauge fixing, for instance along the lines of [31, 30] or by means of orthogonal polynomials as in [5].

An example of a positive function: $f(x) = (1+a x^2) \Phi(b x)$ with $\Phi$ a bump function and $ a = 1/900, b = 1/100 $.

The divided difference $f'[\lambda_k,\lambda_l]$ for this function $f$.

The inverse gauge propagator $f'[\lambda_k,\lambda_l]$ for the $N=61$ smallest eigenvalues of the Dirac operator on the circle (i.e. $ \lambda_k,\lambda_l = -30, -29, \ldots, 30$).

In any case, we are now in a position to consider higher-loop contributions to the spectral action, and, in particular, all one-particle irreducible $n$-point Feynman graphs. Their (possibly divergent) amplitudes form the starting point of the renormalization process of the spectral action.

§.§ Ward identity for the gauge propagator

In addition to the Ward identity (<ref>) for the fermion propagator, we claim that we also have the following Ward identity for the gauge propagator:

\begin{equation}
\label{eq:ward-gauge}
\parbox{2.1cm}{
\hspace{.5cm}
\begin{fmfgraph*}(25,20)
\fmfset{arrow_len}{2mm}
\fmfleft{l1,l2}
\fmfright{r1,r2}
%      \fmf{phantom}{l,v1}
\fmf{fermion,left=.2,tension=.5}{l2,v2,l1}
\fmf{photon}{v2,v3}
\fmf{fermion,left=.2,tension=.5}{r1,v3,r2}
%            \fmf{phantom}{v4,r}
\fmffreeze
\fmftop{b}
\fmf{dashes}{v3,b}
\fmflabel{$a$}{b}
\end{fmfgraph*}}
- \parbox{2.2cm}{
\begin{fmfgraph*}(25,20)
\fmfset{arrow_len}{2mm}
\fmfleft{l1,l2}
\fmfright{r1,r2}
%      \fmf{phantom}{l,v1}
\fmf{fermion,left=.2,tension=.5}{l2,v2,l1}
\fmf{photon}{v2,v3}
\fmf{fermion,left=.2,tension=.5}{r1,v3,r2}
%            \fmf{phantom}{v4,r}
\fmffreeze
\fmftop{b}
\fmf{dashes}{v2,b}
\fmflabel{$a$}{b}
\end{fmfgraph*}}
%  \vspace{.7cm}
\begin{fmfgraph*}(25,20)
\fmfset{arrow_len}{2mm}
\fmfleft{l1,l2}
\fmfright{r1,r2}
%      \fmf{phantom}{l,v1}
\fmf{fermion,left=.2,tension=.5}{l2,v2,l1}
\fmftop{b}
\fmf{photon}{b,bc}
\fmflabel{$\diraccom{a}$}{b}
%            \fmf{phantom}{v4,r}
\fmf{photon}{v2,c1}
\fmf{fermion,left=1,tension=.5}{c2,c1}
\fmf{fermion,left=.5}{c1,bc,c2}
\fmf{photon}{c2,v3}
\fmf{fermion,left=.2,tension=.5}{r1,v3,r2}
\end{fmfgraph*}}
\end{equation}

where every fermion loop adds a minus sign. Indeed, the left-hand side is

\begin{align*}
\wick{ \c Q_{ik} \c Q_{lm} a_{mn} }-   \wick{ a_{im} \c  Q_{mk}  \c Q_{ln}}   & =  G_{ik} \delta_{im} \delta_{kl} a_{mn} - G_{ln} \delta_{mn} \delta_{kl} a_{im} \\
&= ( G_{ik}- G_{nk} )\delta_{kl} a_{in}
\end{align*}

while for the right-hand side we use the defining property of the divided differences to find:

\begin{align*}
-   \wick{ \c Q_{ik} \c Q_{rp}} a_{pq} (\lambda_p -\lambda_q)&\wick{ \c Q_{qr} \c  Q_{ln}} f'[\lambda_p, \lambda_q, \lambda_r]\\
&= - G_{ik} \delta_{ip} \delta_{kr} G_{qr} \delta_{qn} \delta_{rl} a_{pq} (\lambda_p -\lambda_q)f'[\lambda_p, \lambda_q, \lambda_r]\\
& =  G_{ik} G_{nk} \left( f'[\lambda_k, \lambda_n] -  f'[\lambda_i, \lambda_k] \right)  \delta_{kl} a_{in}.
\end{align*}

The two expressions coincide because of the very fact that the free propagator is the inverse of the divided difference.

§.§ Two-point functions at one-loop

The two-point graphs at one-loop
are given in Table <ref>. The external fields $V_1,V_2$ should be assigned to the external legs in all different cyclical manners.

The amplitude for the first graph is given by

\begin{align}
\label{eq:2ptA}
\scalebox{.7}{\parbox{2.5cm}{
\begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
\fmfleft{l}
\fmftop{r1}
\fmfbottom{r2}
\fmf{photon,tension=.4,label=$V_1$}{l,v1}
\fmf{fermion,left=.5,tension=.2}{v1,v2,v3,v1}
\fmf{phantom,tension=.1}{v2,r1}
\fmf{phantom,tension=.1}{v3,r2}
\fmfright{r}
\fmf{photon,tension=.3,label=$V_2$}{r,w1}
\fmf{fermion,left=.5,tension=.2}{w1,w3,w2,w1}
\fmf{phantom,tension=.1}{w2,r1}
\fmf{phantom,tension=.1}{w3,r2}
\fmf{phantom,tension=.05}{v2,w2}
\fmf{phantom,tension=.05}{v3,w3}
\fmfset{wiggly_len}{2mm}
\fmf{photon,tension=.1}{v2,w3}
\fmf{photon,tension=.1}{v3,w2}
\end{fmfgraph*}}} &= \sum_{\begin{smallmatrix} i,j,k \\ l,m,n\end{smallmatrix}} (V_1)_{ij} \wick{ \c1 Q_{jk} \c2 Q_{ki} (V_2)_{lm} \c1 Q_{mn} \c2 Q_{nl} }  f'[\lambda_i,\lambda_j,\lambda_k]  f'[\lambda_l,\lambda_m,\lambda_n]  \nonumber \\
&= \sum_{i,k} (V_1)_{ii} (V_2)_{kk} G_{ik}^2 f'[\lambda_i, \lambda_i,\lambda_k]f'[\lambda_i, \lambda_k,\lambda_k] .
% & \quad = \sum_{i,k} (V_1)_{ii} (V_2)_{kk} \left( \frac{ f''(\lambda_i)- f''(\lambda_k) } { f'(\lambda_i)- f'(\lambda_k) } \right)^2
\end{align}

In particular, there is no running loop index in this expression and so this diagram remains finite even when the size $N$ of the matrices is sent to $\infty$. We conclude that the amplitude of this graph is not relevant for renormalization purposes.

We then turn to the second graph in Table <ref>, and compute

\begin{align}
\scalebox{.7}{ \parbox{2.5cm}{
\begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
\fmfset{wiggly_len}{2mm}
\fmfleft{l}
\fmftop{r1}
\fmfbottom{r2}
\fmf{photon,tension=.4,label=$V_1$}{l,v1}
\fmf{fermion,left=.5,tension=.2}{v1,v2,v3,v1}
\fmf{phantom,tension=.1}{v2,r1}
\fmf{phantom,tension=.1}{v3,r2}
\fmfright{r}
\fmf{photon,tension=.3,label=$V_2$}{r,w1}
\fmf{fermion,left=.5,tension=.2}{w1,w3,w2,w1}
\fmf{phantom,tension=.1}{w2,r1}
\fmf{phantom,tension=.1}{w3,r2}
\fmf{photon,tension=.2}{v2,w2}
\fmf{photon,tension=.2}{v3,w3}
\end{fmfgraph*}}}&= \sum_{\begin{smallmatrix} i,j,k \\ l,m,n\end{smallmatrix}} (V_1)_{ij} \wick{ \c1 Q_{jk} \c2 Q_{ki} (V_2)_{lm} \c2 Q_{mn} \c1 Q_{nl} } f'[\lambda_i,\lambda_j,\lambda_k]  f'[\lambda_l,\lambda_m,\lambda_n]   \nonumber \\
&= \sum_{i,j,k} (V_1)_{ij} (V_2)_{ji} G_{ik}G_{kj} f'[\lambda_i, \lambda_j,\lambda_k]^2.
\label{eq:vertexcontr}
\end{align}

We find that this amplitude has a potential divergence in the limit that $N \to \infty$ (see Figure <ref> for the behaviour of the summands). As such it should be subtracted from the effective action in order to render the theory finite after removal of the regulator.

Summands in (<ref>) for $(\lambda_i,\lambda_j) = (0,0)$

Summands in (<ref>) for $(\lambda_i,\lambda_j) = (10,0)$

Summands in (<ref>) for $(\lambda_i,\lambda_j) = (0,0)$

Summands in (<ref>) for $(\lambda_i,\lambda_j) = (10,0)$

The behaviour of the summands (indexed by $\lambda_k$ running from $-30$ to $30$) for the vertex contribution  in (<ref>) and (<ref>)
for the Dirac operator on the circle and function $f$ as in Figure <ref>.

For the final diagram with two external lines we compute its amplitude to be:

\begin{align}
\scalebox{.7}{  \parbox{3cm}{
\begin{fmfgraph*}(30,30)
\fmfset{arrow_len}{2mm}
\fmfleft{l1}
\fmfright{l2}
\fmftop{r1,r2}
\fmf{photon,tension=1,label=$V_1$}{l1,v1}
\fmf{photon,tension=1,label=$V_2$}{l2,v2}
\fmf{fermion,left=.5,tension=2}{v1,v4,v3,v2,v1}
\fmf{phantom,tension=1}{v4,r1}
\fmf{phantom,tension=1}{v3,r2}
\fmf{photon,left,tension=.05}{v4,v3}
\end{fmfgraph*}}}&= \sum_{ i,j,k,l} (V_1)_{ij} \wick{ \c Q_{jk} \c Q_{kl} } (V_2)_{li}\label{eq:vertexcontr2}
f'[\lambda_i,\lambda_j,\lambda_k,\lambda_l] \\ \nonumber
&= \sum_{i,j,k} (V_1)_{ij} (V_2)_{ji} G_{jk} f'[\lambda_i, \lambda_j,\lambda_j,\lambda_k].\nonumber
\end{align}

Again, this graph amplitude is potentially divergent in the limit $N \to \infty$ and should thus be subtracted. The same applies to the same graph but with $V_1$ and $V_2$ exchanged.

.85 3cm

.85 3cm

The two-point graphs at one-loop.

.85 2.5cm

The relevant three-point graphs at one-loop.

§.§ One-loop counterterms to the spectral action

The computations of the graph amplitudes in the previous section show that the second two graphs in Table <ref> are the relevant ones to consider as counterterms for the spectral action. However, since the spectral action is in particular a gauge theory, it is crucial that such counterterms are of the same form as the terms appearing in the spectral action.

As may be expected, a crucial role will be played by so-called quantum Ward identities. They form the analogue of (<ref>) for the divergent component of the 1PI $n$-point functions at one loop. Let us denote by $\llangle V_1, \ldots, V_n \rrangle^{1L}$ all one-loop $n$-point graphs whose amplitudes involve a sum over a loop index. The skeletons for such graphs are depicted in Table <ref>, for which all external lines are written outside the graph diagram, and labelled in cyclical order. Indeed, if an external line would be in the interior of the diagram, it is surrounded by the loop in the diagram, and will thus prevent the loop index from running (as in Equation <ref>).

The quantum Ward identities are now given by

\begin{multline*}
\dbr{V_1,\ldots,aV_j,\ldots,V_n}^{1L}-\dbr{V_1,\ldots,V_{j-1}a,\ldots,V_n}^{1L}\\
\end{multline*}

It is this identity, in combination with cyclicity of the bracket $\llangle V_1, \ldots, V_n \rrangle = \llangle V_n , V_1, \ldots, V_{n-1} \rrangle$, which allows us to follow line-by-line the derivation of the Chern–Simons and Yang–Mills terms in the previous section (cf. [40]). We thus arrive at our main conclusion which is that

the divergent part of the one-loop quantum effective spectral action  can be expanded as

%   \Gamma_D[V]_{\infty}^{1L}
\sum_n \frac 1n \dbr{V,\ldots,V}_{\infty}^{1L}
= \sum_{k=1}^\infty \left( \int_{\widetilde \psi_{2k-1}}  \!\!\!\!\! \cs_{2k-1} (A) +\frac 1 {2k} \int_{\widetilde \phi_{2k}} \!\!\!\!\! F^{k} \right).

Here $\widetilde \phi$ and $\widetilde \psi$ are the analogues of $\phi$ and $\psi$ as defined in (<ref>) and (<ref>) but now using the double bracket.

We conclude that the passage to the one-loop renormalized spectral action can be realized by a transformation in the space of noncommutative integrals, sending $\phi \mapsto \phi - \widetilde \phi$ and $\psi \mapsto \psi- \widetilde \psi$, thus rendering the theory (one-loop) renormalizable as a gauge theory.

Before addressing the general case of $n$-point vertex contributions, we will present a diagrammatic proof of the quantum Ward identity for divergent one-loop two-point functions.

We first consider the contribution from the second diagram in Table <ref> to the term $\llangle a V_1, V_2\rrangle - \llangle V_1, V_2 a \rrangle $ in the quantum Ward identity:

\begin{align*}
&  \scalebox{.7}{ \parbox{3cm}{
\begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
\fmfleft{lp,l}
\fmftop{r1}
\fmfbottom{r2}
\fmf{photon,tension=.3,label=$V_1$}{l,v1}
\fmf{dashes,tension=.3,label=$a$}{v1,lp}
\fmf{fermion,left=.5,tension=.2}{v1,v2,v3,v1}
\fmf{phantom,tension=.1}{v2,r1}
\fmf{phantom,tension=.1}{v3,r2}
\fmfright{r}
\fmf{photon,tension=.3,label=$V_2$}{r,w1}
\fmf{fermion,left=.5,tension=.2}{w1,w3,w2,w1}
\fmf{phantom,tension=.1}{w2,r1}
\fmf{phantom,tension=.1}{w3,r2}
\fmf{photon,tension=.1}{v2,w2}
\fmf{photon,tension=.1}{v3,w3}
\end{fmfgraph*}}}
\scalebox{.7}{ \parbox{3cm}{
\begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
\fmfleft{l}
\fmftop{r1}
\fmfbottom{r2}
\fmf{photon,tension=.3,label=$V_1$}{l,v1}
\fmf{fermion,left=.5,tension=.2}{v1,v2,v3,v1}
\fmf{phantom,tension=.1}{v2,r1}
\fmf{phantom,tension=.1}{v3,r2}
\fmfright{r,rp}
\fmf{photon,tension=.3,label=$V_2$}{w1,rp}
\fmf{dashes,tension=.3,label=$a$}{w1,r}
\fmf{fermion,left=.5,tension=.2}{w1,w3,w2,w1}
\fmf{phantom,tension=.1}{w2,r1}
\fmf{phantom,tension=.1}{w3,r2}
\fmf{photon,tension=.1}{v2,w2}
\fmf{photon,tension=.1}{v3,w3}
\end{fmfgraph*}}}
\\
&\quad = \quad \scalebox{.7}{  \parbox{3cm}{
\blue{    \begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
\fmfleft{l1,l,l2}
\fmftop{r1}
\fmfbottom{r2}
\fmf{phantom,tension=.05}{l,v1}
\fmf{phantom,tension=.05}{l,v1p}
\fmf{fermion,foreground=blue,left=.5,tension=.2}{v1,v1p,v2,v3,v1}
\fmf{phantom,tension=.1}{v2,r1}
\fmf{phantom,tension=.1}{v3,r2}
\fmfright{r}
\fmf{photon,foreground=blue,tension=.3,label=$V_2$}{r,w1}
\fmf{fermion,foreground=blue,left=.5,tension=.2}{w1,w3,w2,w1}
\fmf{phantom,tension=.1}{w2,r1}
\fmf{phantom,tension=.1}{w3,r2}
\fmf{photon,foreground=blue,tension=.1}{v2,w2}
\fmf{photon,foreground=blue,tension=.1}{v3,w3}
\fmf{photon,foreground=blue,tension=.05,label=$\diraccom{a}$}{v1,l1}
\fmf{photon,foreground=blue,tension=.05,label=$V_1$}{l2,v1p}
\end{fmfgraph*}}}}
+ \quad   \scalebox{.7}{  \parbox{3cm}{
\blue{\begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
\fmfright{l2,l,l1}
\fmfbottom{r1}
\fmftop{r2}
\fmf{phantom,tension=.05}{l,v1}
\fmf{phantom,tension=.05}{l,v1p}
\fmf{fermion,foreground=blue,left=.5,tension=.2}{v1,v1p,v2,v3,v1}
\fmf{phantom,tension=.1}{v2,r1}
\fmf{phantom,tension=.1}{v3,r2}
\fmfleft{r}
\fmf{photon,foreground=blue,tension=.3,label=$V_1$}{r,w1}
\fmf{fermion,foreground=blue,left=.5,tension=.2}{w1,w3,w2,w1}
\fmf{phantom,tension=.1}{w2,r1}
\fmf{phantom,tension=.1}{w3,r2}
\fmf{photon,foreground=blue,tension=.1}{v2,w2}
\fmf{photon,foreground=blue,tension=.1}{v3,w3}
\fmf{photon,foreground=blue,tension=.05,label=$V_1$}{v1,l1}
\fmf{photon,foreground=blue,tension=.05,label=$\diraccom{a}$}{l2,v1p}
\end{fmfgraph*}}}}
%  \\[4mm]   &
\scalebox{.7}{ \parbox{2.2cm}{
\green{    \begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
\fmftop{l,t,r}
\fmf{photon,foreground=green,tension=.7,label=$V_1$}{l,v1}
\fmf{fermion,foreground=green,left=.5,tension=.3}{v1,v2p,v2,v1}
\fmf{phantom,tension=.2}{v2p,t}
\fmfbottom{b}
\fmf{photon,foreground=green,tension=.7,label=$\diraccom{a}$}{b,bc}
\fmf{photon,foreground=green,tension=.7,label=$V_2$}{v4,r}
\fmf{fermion,foreground=green,left=.5,tension=.3}{c2,bc,c1,c2}
\fmf{photon,foreground=green,tension=.4}{v2,c1}
\fmf{photon,foreground=green,tension=.4}{c2,v3}
\fmf{fermion,foreground=green,left=.5,tension=.3}{v3,v4p,v4,v3}
\fmf{phantom,tension=.2}{v4p,t}
\fmf{photon,foreground=green,tension=.4}{v2p,v4p}
\end{fmfgraph*}}}}
\end{align*}

For the third two-point diagram in Table <ref> there are two possible assignments of the external fields, so that their contribution to $\llangle a V_1, V_2\rrangle - \llangle V_1, V_2 a \rrangle $ is

\begin{align*}
\scalebox{.7}{ \parbox{2.2cm}{
\begin{fmfgraph*}(20,20)
\fmfset{arrow_len}{2mm}
\fmfleft{l1}
\fmfright{l2}
\fmftop{r1,r2}
\fmf{photon,tension=1,label=$V_1$}{l1,v1}
\fmf{photon,tension=1,label=$V_2$}{l2,v2}
\fmf{fermion,left=.5,tension=2}{v1,v4,v3,v2,v1}
\fmf{phantom,tension=1}{v4,r1}
\fmf{phantom,tension=1}{v3,r2}
\fmf{photon,left,tension=.05}{v4,v3}
\fmffreeze
\fmfbottom{b}
\fmf{dashes}{v1,b}
\fmflabel{$a$}{b}
\end{fmfgraph*}}}  -
\quad
\scalebox{.7}{  \parbox{2.5cm}{
\begin{fmfgraph*}(20,20)
\fmfset{arrow_len}{2mm}
\fmfleft{l1}
\fmfright{l2}
\fmftop{r1,r2}
\fmf{photon,tension=1,label=$V_1$}{l1,v1}
\fmf{photon,tension=1,label=$V_2$}{l2,v2}
\fmf{fermion,left=.5,tension=2}{v1,v4,v3,v2,v1}
\fmf{phantom,tension=1}{v4,r1}
\fmf{phantom,tension=1}{v3,r2}
\fmf{photon,left,tension=.05}{v4,v3}
\fmffreeze
\fmfbottom{b}
\fmf{dashes}{v2,b}
\fmflabel{$a$}{b}
\end{fmfgraph*}}}
=  \scalebox{.7}{ \parbox{3cm}{
\vspace{-1.5cm}
\red{ \begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
%% \fmfleft{l1}
%% \fmfright{l2}
\fmftop{r1,r2}
\fmfbottom{l1,b,l2}
\fmf{photon,foreground=red,tension=1,label=$V_1$}{l1,v1}
\fmf{photon,foreground=red,tension=1,label=$V_2$}{l2,v2}
\fmf{fermion,foreground=red,left=.3,tension=3}{v1,v4,v3,v2,bc,v1}
\fmf{phantom,tension=1}{v4,r1}
\fmf{phantom,tension=1}{v3,r2}
\fmf{photon,foreground=red,left,tension=.01}{v4,v3}
\fmf{photon,foreground=red}{b,bc}
\fmflabel{$\diraccom{a}$}{b}
\end{fmfgraph*}}}}
\end{align*}

\begin{align*}
&   \scalebox{.7}{ \parbox{3cm}{
\begin{fmfgraph*}(20,20)
\fmfset{arrow_len}{2mm}
\fmfleft{l1}
\fmfright{l2}
\fmftop{r1,r2}
\fmf{photon,tension=1,label=$V_2$}{v1,l1}
\fmf{photon,tension=1,label=$V_1$}{v2,l2}
\fmf{fermion,left=.5,tension=2}{v1,v4,v3,v2,v1}
\fmf{phantom,tension=1}{v4,r1}
\fmf{phantom,tension=1}{v3,r2}
\fmf{photon,left,tension=.05}{v4,v3}
\fmffreeze
\fmf{dashes}{v2,r2}
\fmflabel{$a$}{r2}
\end{fmfgraph*}} }-
\quad    \scalebox{.7}{ \parbox{2.2cm}{
\begin{fmfgraph*}(20,20)
\fmfset{arrow_len}{2mm}
\fmfleft{l1}
\fmfright{l2}
\fmftop{r1,r2}
\fmf{photon,tension=1,label=$V_2$}{v1,l1}
\fmf{photon,tension=1,label=$V_1$}{v2,l2}
\fmf{fermion,left=.5,tension=2}{v1,v4,v3,v2,v1}
\fmf{phantom,tension=1}{v4,r1}
\fmf{phantom,tension=1}{v3,r2}
\fmf{photon,left,tension=.05}{v4,v3}
\fmffreeze
\fmf{dashes}{v1,r1}
\fmflabel{$a$}{r1}
\end{fmfgraph*}}}\\
& \quad =  \scalebox{.7}{ \parbox{3cm}{
\vspace{-1.5cm}
\red{   \begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
%% \fmfleft{l1}
%% \fmfright{l2}
\fmftop{r1,r2}
\fmfbottom{l1,b,l2}
\fmf{photon,foreground=red,tension=1,label=$V_2$}{l1,v1}
\fmf{photon,foreground=red,tension=1,label=$\diraccom{a}$}{l2,v2}
\fmf{fermion,foreground=red,left=.3,tension=3}{v1,v4,v3,v2,bc,v1}
\fmf{phantom,tension=1}{v4,r1}
\fmf{phantom,tension=1}{v3,r2}
\fmf{photon,foreground=red,left,tension=.01}{v4,v3}
\fmf{photon,foreground=red,label=$V_1$}{bc,b}
\end{fmfgraph*}}}}
+ \quad \scalebox{.7}{  \parbox{2.7cm}{
\vspace{-1.5cm}
\red{   \begin{fmfgraph*}(25,25)
\fmfset{arrow_len}{2mm}
%% \fmfleft{l1}
%% \fmfright{l2}
\fmftop{r1,r2}
\fmfbottom{l1,b,l2}
\fmf{photon,foreground=red,tension=1,label=$\diraccom{a}$}{l1,v1}
\fmf{photon,foreground=red,tension=1,label=$V_1$}{l2,v2}
\fmf{fermion,foreground=red,left=.3,tension=3}{v1,v4,v3,v2,bc,v1}
\fmf{phantom,tension=1}{v4,r1}
\fmf{phantom,tension=1}{v3,r2}
\fmf{photon,foreground=red,left,tension=.01}{v4,v3}
\fmf{photon,foreground=red,label=$V_2$}{bc,b}
\end{fmfgraph*}}}}
%\\ & \qquad \qquad
+    \scalebox{.7}{ \parbox{3cm}{
\blue{ \begin{fmfgraph*}(35,20)
\fmfset{arrow_len}{2mm}
\fmfleft{l1,l,l2}
\fmftop{r1}
\fmfbottom{r2}
\fmf{phantom,tension=.05}{l,v1}
\fmf{phantom,tension=.05}{l,v1p}
\fmf{fermion,foreground=blue,left=.5,tension=.2}{v1,v1p,v2,v3,v1}
\fmf{phantom,tension=.1}{v2,r1}
\fmf{phantom,tension=.1}{v3,r2}
\fmfright{r}
\fmf{photon,foreground=blue,tension=.3,label=$\diraccom{a}$}{r,w1}
\fmf{fermion,foreground=blue,left=.5,tension=.18}{w1,w3,w2,w1}
\fmf{phantom,tension=.1}{w2,r1}
\fmf{phantom,tension=.1}{w3,r2}
\fmf{photon,foreground=blue,tension=.2}{v2,w2}
\fmf{photon,foreground=blue,tension=.2}{v3,w3}
\fmf{photon,foreground=blue,tension=.05,label=$V_1$}{l1,v1}
\fmf{photon,foreground=blue,tension=.05,label=$V_2$}{v1p,l2}
\end{fmfgraph*}}}}
\end{align*}

We have coloured the Feynman graphs on the right-hand side of the quantum Ward identity according to their topology, i.e. as they appear in Table <ref>. One then readily sees that the graphs conspire to yield all cyclic permutations of $\diraccom{a},V_1,V_2$ as external fields on all planar one-loop graphs with three external legs.

dotssubpath (length(p)/6,5*length(p)/6) of p

dotssubpath (length(p1)/6,5*length(p1)/6) of p1
dotssubpath (length(p2)/6,5*length(p2)/6) of p2

dotssubpath (length(p1)/6,5*length(p1)/6) of p1
dotssubpath (length(p2)/6,5*length(p2)/6) of p2
dotssubpath (length(p3)/6,5*length(p1)/6) of p3

dotssubpath (length(p1)/6,5*length(p1)/6) of p1
dotssubpath (length(p2)/6,5*length(p2)/6) of p2
dotssubpath (length(p3)/6,5*length(p3)/6) of p3
dotssubpath (length(p4)/6,5*length(p4)/6) of p4

Skeletons for divergent one-loop $n$-point functions with increasing number of vertices. The fermion loops that define the vertices are all oriented as clockwise.

This argument extends to all potentially divergent one-loop $n$-point functions $\llangle V_1, \ldots, V_n \rrangle^{1L}$ as follows. Recall that all such divergent one-loop diagrams have skeletons as depicted in Table <ref>, with the external lines labelled cyclically from $1$ to $n$.

The decoration of the external legs of our graphs with the external fields $V_1, \ldots, V_n$ then proceeds according to this labelling $1, \ldots, n$ and, upon summing over all such decorated graphs $G$, we get

\llangle V_1, \ldots, V_n \rrangle^{1L} = \sum_G %\sum_{\sigma \in C_n} \frac 1 {|\Aut(G)|}
G_{ V_{1},\ldots V_{n} }.

The left-hand side of the quantum Ward identity essentially comes down to connecting external edges to the graphs $G$. We will write $G_i$ for the graph $G$ with an insertion of an external gauge edge at a point $i$ in between $n$ and $1$: this insertion point $i$ can be either an outer fermion line in $G$ (as in (<ref>)) or, if $1$ and $n$ are not attached to the same vertex in $G$, a gauge propagator (as in (<ref>)). We then find

\begin{align*}
\llangle a  V_1, \ldots, V_n \rrangle^{1L} -
\llangle V_1, \ldots, V_n  a \rrangle^{1L} & = \sum_{G,i} %\sum_{i} \frac 1 {|\Aut(G)|}
(G_i)_{ \diraccom{a} , V_1, \ldots ,V_n  },
\end{align*}

where the decoration $\diraccom{a}$ is attached to the external gauge edge inserted at the point $i$ of $G_i$.

It is clear that the sum over $G$ and $i$ yield all decorated $n+1$-point graphs, and, moreover, that any $n+1$-point graph with labels $\diraccom{a} , V_1, \ldots ,V_n$ is obtained in a unique manner from an insertion of an external edge in an $n$-point graph, as described above. We are thus left precisely with $\llangle {\diraccom{a}},  V_1, \ldots, V_n \rrangle^{1L}$ as desired.

§ CONCLUSIONS

In this paper we have analyzed the quantum gauge fluctuations for the spectral action in noncommutative geometry. Using the background field method we have showed one-loop renormalizability of the spectral action, while staying within the same spectral framework.

Naturally, this forms the starting point for more direct applications of noncommutative geometry to particle physics phenomenology. Instead of the spectral action playing the role of a bare action functional, to which subsequent RG-methods are applied, we now have a candidate for a so-called quantum effective spectral action, given by the sum of all 1PI Feynman diagrams and which is supposed to be valid at all energies.
One may then try to extend the derivation of bare physical Lagrangians from the spectral action [9, 39] to the renormalized spectral action, and arrive at a spectral, noncommutative geometric description of particle physics which is also valid and falsifiable at lower energies.

Besides these future steps in the applications to particle physics phenomenology, it is also important to extend the “power-counting” and diagrammatics of the one-loop renormalizability that we presented here to arbitrary loop order. This, and also a more detailed account of the derivation presented in this paper, will be reported elsewhere.
The connection with the proof of renormalizability for noncommutative scalar field theories [28] also deserves further investigation. One of the main differences is that they consider so-called non-local matrix models [27] with a quartic vertex, while instead we have a local matrix model but with vertices of arbitrary valence.

§.§ Acknowledgements

We thank Steven Lord, Fedor Sukochev and the other members in their research group for fruitful discussions during a visit in August 2019, where the basis for the current paper was laid. We also thank Ali Chamseddine and Alain Connes for useful comments.

Research supported by NWO Physics Projectruimte (680-91-101).

[1]
N. Alkofer, F. Saueressig, and O. Zanusso.
Spectral dimensions from the spectral action.
Phys. Rev. D91 (2015)  025025.

[2]
S. Azarfar and M. Khalkhali.
Random finite noncommutative geometries and topological recursion,

[3]
J. W. Barrett and L. Glaser.
Monte Carlo simulations of random non-commutative geometries.
J. Phys. A49 (2016)  245001.

[4]
W. Beenakker, T. van den Broek, and W. D. van Suijlekom.
Supersymmetry and noncommutative geometry, volume 9 of SpringerBriefs in Mathematical Physics.
Springer, Cham, 2016.

[5]
D. Bessis, C. Itzykson, and J. B. Zuber.
Quantum field theory techniques in graphical enumeration.
Adv. Appl. Math. 1 (1980)  109–157.

[6]
A. Bochniak and A. Sitarz.
Spectral geometry for the standard model without fermion doubling.
Phys. Rev. D 101 (2020)  075038.

[7]
L. Boyle and S. Farnsworth.
Non-Commutative Geometry, Non-Associative Geometry and the Standard
Model of Particle Physics.
New J. Phys. 16 (2014)  123027.

[8]
T. broekvan den Broek and W. D. van Suijlekom.
Supersymmetric QCD from noncommutative geometry.
Phys. Lett. B699 (2011)  119–122.

[9]
A. H. Chamseddine and A. Connes.
Universal formula for noncommutative geometry actions: Unifications
of gravity and the Standard Model.
Phys. Rev. Lett. 77 (1996)  4868–4871.

[10]
A. H. Chamseddine, A. Connes, and M. Marcolli.
Gravity and the Standard Model with neutrino mixing.
Adv. Theor. Math. Phys. 11 (2007)  991–1089.

[11]
A. H. Chamseddine, A. Connes, and V. Mukhanov.
Quanta of geometry: Noncommutative aspects.
Phys. Rev. Lett. 114 (2015)  091302.

[12]
A. H. Chamseddine, A. Connes, and W. D. van Suijlekom.
Beyond the spectral Standard Model: Emergence of Pati-Salam
JHEP 1311 (2013)  132.

[13]
A. H. Chamseddine, A. Connes, and W. D. van Suijlekom.
Grand unification in the spectral Pati-Salam model.
JHEP 11 (2015)  011.

[14]
A. H. Chamseddine, J. Iliopoulos, and W. D. van Suijlekom.
Spectral action in matrix form.
Eur. Phys. J. C 80 (2020)  1045.

[15]
A. Connes.
Noncommutative Geometry.
Academic Press, San Diego, 1994.

[16]
A. Connes.
Gravity coupled with matter and the foundation of non-commutative
Commun. Math. Phys. 182 (1996)  155–176.

[17]
A. Connes and M. Marcolli.
Noncommutative Geometry, Quantum Fields and Motives.
AMS, Providence, 2008.

[18]
A. Connes and A. H. Chamseddine.
Inner fluctuations of the spectral action.
J. Geom. Phys. 57 (2006)  1–21.

[19]
L. Dabrowski, F. D'Andrea, and A. Sitarz.
The Standard Model in noncommutative geometry: fundamental
fermions as internal forms.
Lett. Math. Phys. 108 (2018)  1323–1340.

[20]
L. Dabrowski and A. Sitarz.
Fermion masses, mass-mixing and the almost commutative geometry of
the Standard Model.
JHEP 02 (2019)  068.

[21]
A. Devastato, F. Lizzi, and P. Martinetti.
Higgs mass in noncommutative geometry.
Fortschr. Phys. 62 (2014)  863–868.

[22]
A. Devastato, F. Lizzi, and P. Martinetti.
Grand symmetry, spectral action, and the Higgs mass.
JHEP 1401 (2014)  042.

[23]
A. Devastato and P. Martinetti.
Twisted spectral triple for the Standard Model and spontaneous
breaking of the Grand Symmetry.
Math. Phys. Anal. Geom. 20 (2017)  2.

[24]
L. Glaser and A. Stern.
Understanding truncated non-commutative geometries through computer
J. Math. Phys. 61 (2020)  033507.

[25]
L. Glaser and A. B. Stern.
Reconstructing manifolds from truncations of spectral triples.
J. Geom. Phys. 159 (2021)  Paper No. 103921, 17.

[26]
J. Gomis and S. Weinberg.
Are nonrenormalizable gauge theories renormalizable?
Nucl. Phys. B469 (1996)  473–487.

[27]
H. Grosse and R. Wulkenhaar.
Power-counting theorem for non-local matrix models and
Comm. Math. Phys. 254 (2005)  91–127.

[28]
H. Grosse and R. Wulkenhaar.
Renormalization of $\phi^4$ theory on noncommutative $\mathbb R^4$
in the matrix base.
Commun. Math. Phys. 256 (2005)  305–374.

[29]
B. Iochum, C. Levy, and D. Vassilevich.
Spectral action beyond the weak-field approximation.
Commun. Math. Phys. 316 (2012)  595–613.

[30]
R. A. Iseppi.
The BV formalism: theory and application to a matrix model.
Rev. Math. Phys. 31 (2019)  1950035, 24.

[31]
R. A. Iseppi and W. D. van Suijlekom.
Noncommutative geometry and the BV formalism: application to a
matrix model.
J. Geom. Phys. 120 (2017)  129–141.

[32]
M. Khalkhali and N. Pagliaroli.
Phase transition in random noncommutative geometries.
Journal of Physics A: Mathematical and Theoretical 54 (2020)

[33]
M. A. Kurkov, F. Lizzi, and D. Vassilevich.
High energy bosons do not propagate.
Phys. Lett. B 731 (2014)  311–315.

[34]
M. Nakahara.
Geometry, Topology and Physics.
IOP Publishing, 1990.

[35]
A. Sitarz.
Towards the signs of new physics through the spectral action.
Int. J. Geom. Meth. Mod. Phys. 17 (2020)  2040008.

[36]
A. Skripka.
Asymptotic expansions for trace functionals.
J. Funct. Anal. 266 (2014)  2845–2866.

[37]
W. D. suijlekomvan Suijlekom.
Perturbations and operator trace functions.
J. Funct. Anal. 260 (2011)  2483–2496.

[38]
W. D. suijlekomvan Suijlekom.
Renormalization of the spectral action for the Yang-Mills system.
JHEP 1103 (2011)  146.

[39]
W. D. suijlekomvan Suijlekom.
Noncommutative Geometry and Particle Physics.
Springer, 2015.

[40]
T. D. H. van Nuland and W. D. van Suijlekom.
Cyclic cocycles in the spectral action.
J. Noncommut. Geom., to appear (arXiv:2104.09899).


# Modeling material transport regulation and traffic jam in neurons using PDE-
constrained optimization

Angran Li
Department of Mechanical Engineering
Carnegie Mellon University
Pittsburgh, PA 15213
<EMAIL_ADDRESS>
& Yongjie Jessica Zhang*
Department of Mechanical Engineering &
Department of Biomedical Engineering
Carnegie Mellon University
Pittsburgh, PA 15213
<EMAIL_ADDRESS>

###### Abstract

The intracellular transport process plays an important role in delivering
essential materials throughout branched geometries of neurons for their
survival and function. Many neurodegenerative diseases have been associated
with the disruption of transport. Therefore, it is essential to study how
neurons control the transport process to localize materials to necessary
locations. Here, we develop a novel optimization model to simulate the traffic
regulation mechanism of material transport in complex geometries of neurons.
The transport is controlled to avoid traffic jam of materials by minimizing a
pre-defined objective function. The optimization subjects to a set of partial
differential equation (PDE) constraints that describe the material transport
process based on a macroscopic molecular-motor-assisted transport model of
intracellular particles. The proposed PDE-constrained optimization model is
solved in complex tree structures by using isogeometric analysis (IGA).
Different simulation parameters are used to introduce traffic jams and study
how neurons handle the transport issue. Specifically, we successfully model
and explain the traffic jam caused by reduced number of microtubules (MTs) and
MT swirls. In summary, our model effectively simulates the material transport
process in healthy neurons and also explain the formation of a traffic jam in
abnormal neurons. Our results demonstrate that both geometry and MT structure
play important roles in achieving an optimal transport process in neuron.

_Keywords_ Neuron transport $\cdot$ Traffic jam $\cdot$ Microtubule swirls
$\cdot$ PDE-constrained optimization $\cdot$ Isogeometric analysis

## Introduction

The neuron exhibits a highly polarized structure that typically consists of a
single long axon and multiple dendrites which are both extended from its cell
body. Since most of the materials necessary for the neuron are synthesized in
the cell body, they need to experience long-distance transport in axons or
dendrites to reach their effective location [1, 2]. The intracellular material
transport is therefore especially crucial to ensure necessary materials are
delivered to the right locations for the development, function, and survival
of neuron cells. The disruption of intracellular transport can lead to the
abnormal accumulations of certain cellular material and extreme swelling of
the axon, which have been observed in many neurological and neurodegenerative
diseases such as Huntington’s, Parkinson’s, and Alzheimer’s disease [3, 4, 5,
6, 7]. Therefore, it is essential to study and understand mechanisms of the
transport function and dysfunction.

Recent studies have shown that the neuron is critically dependent on molecular
motors to transport various materials along the longitudinal cytoskeletal
structure like microtubules (MTs) [8, 9, 10]. MTs are long and polarized
polymers with biophysically distinct plus and minus ends [11, 12, 13]. The
polarity of MTs can decide the preferred direction in which an individual
molecular motor moves. For instance, molecular motors from the kinesin and
dynein superfamilies have been identified to convey materials along MTs
towards their plus and minus ends respectively [14, 15]. Inspired by these
findings, there have been many mathematical models proposed to quantitatively
study the motor-driven transport process and understand the pathology of
neuron diseases. For instance, the partial differential equations (PDEs) of
linear reaction-hyperbolic form have been used to approximate the traveling
waves of a single moving species [16]. This model was further extended to
account for multiple moving species [17] and their diffusion [18, 19]. Based
on PDE-based transport, stochastic models have also been developed for both
axonal transport [20, 21] and dendritic transport [22, 23]. In addition,
several mathematical models were developed to simulate material transport in
unhealthy neurons. Xue et al. presented a stochastic model to explain the
segregation of MTs and neurofilaments in neurological diseases [24]. Bertsch
et al. proposed to couple Smoluchowski equations and kinetic-type transport
equations to study the onset and progression of Alzheimer’s disease [25].

Though the aforementioned PDE and stochastic models can successfully simulate
and explain certain phenomena during transport, most of these models were
solved only in simple one-dimensional (1D) or 2D domains without considering
the complex neuron geometry. Recent developments in numerical methods allow us
to obtain accurate solution of PDEs in complex geometries. Specifically,
isogeometric analysis (IGA) [26] directly integrates geometric modeling with
numerical simulation and achieves better accuracy and robustness compared to
the conventional finite element method (FEM), making it a perfect tool to
tackle the highly branched neuron geometry. In particular, IGA performs
simulation with different types of splines as basis functions instead of
Lagrange polynomials used in conventional FEM. The same smooth spline basis
functions [27] used for both geometrical modeling and numerical simulation
lead to accurate geometry representation with high-order continuity and
superior numerical accuracy in simulation. Therefore, IGA has been extensively
used in shell analysis [28, 29, 30, 31], cardiovascular modeling [32, 33, 34,
35, 36, 37], neuroscience simulation [38, 39], fluid-structure interaction
[40, 41, 42, 43], as well as industrial application [44, 45]. Truncated
T-splines [46, 47] were developed to support local refinement over
unstructured quadrilateral and hexahedral meshes. Blended B-splines [48] and
Catmull-Clark subdivision basis functions [49] were investigated to enable
improved or even optimal convergence rates for IGA. With the advances in IGA,
we developed an IGA-based simulation platform to accurately reconstruct
complex neuron geometries and solved a 3D motor-assisted transport model
within them [38]. We also developed a deep learning framework based on the IGA
simulation platform to predict the material transport process in complex
neurite networks [50]. The results from our IGA solver showed how the complex
neuron geometry affects the spatiotemporal material distribution at neurite
junctions and within different branches. However, the motor-assisted model
only provides a simplified model of the actual transport process but ignores
the active regulation from neuron itself.

To model the active regulation from neurons to control the transport process,
we propose to use PDE-constrained optimization (PDE-CO). PDEs are commonly
used in science and engineering to mathematically represent biological and
physical phenomena. Recent advances in numerical methods and high-performance
computing equip the development of large-scale PDE solvers. As a result, PDE-
CO problems arise in a variety of applications including optimal design [51,
52, 53], optimal control [54, 55, 56], and inverse problem [57, 58]. In
particular, PDE-CO has important biomedical applications in exploiting
valuable information from real medical data. For instance, Hogea et al.
presented a PDE-CO framework for modeling gliomas growth and their mass-effect
on the surrounding brain tissue [59]. Kim et al. proposed a transport-theory-
based PDE-constrained multispectral imaging algorithm to reconstruct the
spatial distribution of chromophores in tissue [60]. Melani utilized the blood
flow data and solved a PDE-CO problem based on fluid-structure interaction to
estimate the compliance of arterial walls in vascular networks [61]. PDE-CO
problems was also used to model tumor growth model by fitting the numerical
solution with real experiment data and estimating unknown parameters in the
model [62, 63].

In this study, we develop a novel IGA-based PDE-CO framework to simulate the
material transport regulation and investigate the formation of traffic jams
and swirl during the transport process in complex neurite structures.
Specifically, we design a new objective function in the PDE-CO model to
simulate the control mechanism to (1) mediate the transport velocity field;
and (2) avoid the traffic jam caused by local material accumulation. The
control strength can be adjusted through two penalty parameters in the
objective function. We can also modify the governing PDEs to study the
formation of traffic jam. In particular, our model can simulate the traffic
jam caused by the reduction of MTs and MT swirls during transport. To further
study the influence of geometry on transport, we simulate material transport
in two neuron tree structures with swelling geometry. In summary, our
simulation reveals that the molecular motors and MT structure play fundamental
roles in controlling the delivery of material by mediating the transport
velocity on MTs. The defective transport on MTs can cause material
accumulation in a local region which may further lead to the degeneration of
neuron cells. Combined with geometry of the neurite network, the motor-
assisted transport on MTs controls the routing of material transport at
junctions of neurite branches and effectively distributes transported
materials throughout the networks. Therefore, our study provides key insights
into how material transport in neurite networks is mediated by MTs and their
complex geometry. Our IGA optimization framework is also transformative and
can be extended to solve other PDE-CO models of cellular processes in complex
neurite networks.

## IGA-based material transport optimization in neurons

Our interest lies in the transport of particles along an axon or dendrite in
neuron cell. In our previous work, we simulated the material transport process
using a macroscopic molecular-motor-assisted transport model without any
transport control [38]. Built upon this transport model, we propose a novel
transport optimization model to further study the transport control mechanism
of neuron and predict the formation of a traffic jam in abnormal neurons. The
proposed optimization problem is described as

minimize
$\displaystyle\mathcal{J}(n_{\pm},\bm{v}_{\pm},f_{\pm})=\frac{1}{2}\int_{0}^{T}\int_{\Omega}(\bm{v}_{\pm}-V_{\pm})^{2}d\Omega
dt+\frac{\alpha}{2}\int_{0}^{T}\int_{\Omega}\left||\nabla
n_{\pm}\right||^{2}d\Omega
dt+\frac{\beta}{2}\int_{0}^{T}\int_{\Omega}f_{\pm}^{2}d\Omega dt,$
$\displaystyle\text{in $\Omega$},$ (1a) subject to
$\displaystyle\dfrac{\partial n_{0}}{\partial
t}=-(k_{+}+k_{-})n_{0}+k^{\prime}_{+}l_{+}n_{+}+k^{\prime}_{-}l_{-}n_{-},$
$\displaystyle\text{in $\Omega$},$ (1b)
$\displaystyle\dfrac{\partial(l_{\pm}n_{\pm})}{\partial
t}+\bm{v}_{\pm}\cdot\nabla(l_{\pm}n_{\pm})=D_{\pm}\nabla^{2}(l_{\pm}n_{\pm})+k_{\pm}n_{0}-k^{\prime}_{\pm}(l_{\pm}n_{\pm}),$
$\displaystyle\text{in $\Omega$},$ (1c)
$\displaystyle\dfrac{\partial\bm{v}_{\pm}}{\partial
t}+\bm{v}_{\pm}\cdot\nabla\bm{v}_{\pm}=-\nabla
n_{\pm}+\nabla\cdot(\mu\nabla\bm{v}_{\pm})+f_{\pm},$ $\displaystyle\text{in
$\Omega$},$ (1d) $\displaystyle n_{0}=n_{i},n_{+}=\lambda_{i}n_{i}$
$\displaystyle\text{at inlet end},$ (1e) $\displaystyle
n_{0}=n_{o},n_{-}=\lambda_{o}n_{o}$ $\displaystyle\text{at outlet end},$ (1f)

where the open set $\Omega\subset R^{d}$ ($d=2$ or 3) represents the
$d$-dimensional internal space of the neuron, $V_{\pm}$ is a predefined
velocity field inside neuron; $n_{0}$ , $n_{+}$ and $n_{-}$ are the spatial
concentrations of free, incoming (relative to the cell body; retrograde), and
outgoing (anterograde) particles, respectively; $D_{\pm}$ is the diffusion
coefficient of incoming and outgoing materials; $\bm{v}_{+}$ and $\bm{v}_{-}$
are velocities of incoming and outgoing particles, respectively; $k_{\pm}$ and
$k^{\prime}_{\pm}$ are rates of MT attachment and detachment of incoming and
outgoing materials, respectively; $l_{\pm}$ represents the density of MTs used
for motor-assisted transport; $f_{\pm}$ represents the control forces that
mediate the material transport; $\mu$ is viscosity of traffic flow;
$\lambda_{i}$, $\lambda_{o}$ represent the degree of loading at inlet and
outlet ends, respectively [18]; and $n_{i}$, $n_{o}$ represent the boundary
value of $n_{0}$ at inlet and outlet ends, respectively. Note that in this
PDE-CO model, $n_{0}$, $n_{\pm}$ and $\bm{v}_{\pm}$ are referred as the “state
variables” while $f_{\pm}$ are referred as the “control variables”. In this
study, we assume the MT system is unipolar that leads to a unidirectional
material transport process and ignore $n_{-}$, $l_{-}$, $\bm{v}_{-}$, $k_{-}$,
$k^{\prime}_{-}$ terms in Eq. 1b-1f. The default values of simulation
parameters are summarized in Table 1.

Herein, we account for active regulation from neuron in the objective function
(Eq. 1a), and we assume the optimal material transportation within neuron can
be achieved by solving the proposed optimization model. The first term in Eq.
1a measures the difference between $\bm{v}_{\pm}$ and the predefined optimal
velocity field $V_{\pm}$. It serves as a velocity control mechanism that
neuron expects to achieve the predefined velocity field $V_{\pm}$ during
transport. The second term measures the cost from concentration gradient
$\nabla n_{\pm}$ within the entire neuron cell. It serves as a traffic jam
control mechanism that the neuron can improve local traffic jam by detecting
and avoiding high concentration gradient in the entire geometry. The value of
parameter $\alpha$ represents to what extent we want to optimize the transport
process and avoid traffic jams. The third term is a regularization that
measures the control forces applied by neuron to mediate the transport. The
value of parameter $\beta$ represents how much the neuron can affect the
transport velocity. To introduce traffic jams in neurons, we modify the
simulation parameters in the governing equations. In this study, we modify the
spatial distribution of $l_{\pm}$ to model the traffic jam caused by abnormal
MTs such as the reduction of MTs and MT swirls during transport.

We employ the “all-at-once” method [64, 65] and IGA to formulate and solve the
optimization model (Eq. 1a) with PDE constraints (Eq. 1b-1f) simultaneously.
We first discretize the objective function to obtain its approximation
$\mathcal{J}_{h}(n_{\pm},\bm{v}_{\pm},f_{\pm})$. We also discretize PDE
constraints (Eq. 1b-1f) to obtain their weak form
$\mathcal{B}_{h}(n_{0},n_{\pm},\bm{v}_{\pm},f_{\pm})=0$. Then, we build a
discrete Lagrangian

$\mathcal{L}_{h}=\mathcal{J}_{h}+p^{T}\mathcal{B}_{h}=0,$ (2)

where $p^{T}$ is the Lagrange multiplier and is also referred to as the
“adjoint variable”. By taking derivatives of the discrete Lagrangian with
respect to state, control, and adjoint variables and setting the resulting
expressions to zero, we obtain the first-order conditions, or Karush-Kuhn-
Tucker (KKT) conditions. The resulting KKT system is then solved using the
GMRES [66] solver implemented in PETSc [67]. In this study, we focus on
solving the proposed optimization model in 2D neuron geometries.

Figure 1: An overview of the material transport control simulation in a
bifurcation geometry. The traffic jam is introduced by reducing MTs in the red
dashed circle region. Color bars unit for velocity field: $\mu m/s$ and
concentration: $mol/\mu m^{3}$.

As shown in Fig. 1, we use a bifurcation example to illustrate the pipeline of
our simulation. We first generate a control mesh and reconstruct the neuron
geometry with Truncated Hierarchical B-splines (THB-spline) by utilizing the
geometry information stored in a SWC file. The SWC file is widely used to
store neuron morphologies including vertices and the associated diameters on
the skeleton of the neuron. We can obtain the SWC files for various real
neuron geometries from the NeuroMorpho database [68]. The raw SWC file needs
to be pre-processed to ensure no duplicated vertices or overlapping skeleton
exist in the geometry. During the geometric modeling of our workflow, we take
the cleaned-up neuron skeleton as input and use the skeleton-based sweeping
method [32] to generate quadrilateral control mesh of the neuron geometry.
Then, we build THB-spline on the quadrilateral mesh [30, 31] for the final
representation of the neuron geometry. Once the spline information for the
geometry is obtained, we run a steady-state Navier-Stokes solver to generate
the pre-defined velocity for the optimization. We then use the default
simulation parameters in Table 1 and modify the spatial distribution of
$l_{+}$ in the red circle regions to introduce traffic jam. Finally, we run
the optimization solver and obtain the velocity field and concentration
distribution. In this paper, we apply the pipeline to various neural
structures with material transport regulation, traffic jam and MT swirl. All
simulations are conducted on the XSEDE (Extreme Science and Engineering
Discovery Environment) supercomputer Bridges at the Pittsburgh Supercomputer
Center [69, 70].

## Results

### Simulation of material transport regulation and traffic jam

We first simulate the normal material transport and the abnormal transport
with traffic jam in a single pipe geometry (Fig. 2). The predefined velocity
field for both cases is computed by solving a steady-state Navier-Stokes
equation and the result is shown in Fig. 2A. The other simulation parameter
settings are summarized in Table 1. The computed velocity field and the
distribution of concentration in the normal transport are shown in Fig. 2B&E.
To model traffic jam caused by the reduction of MTs, the distribution of
$l_{+}$ along the pipe is defined as shown in Fig. 2D. The velocity field and
material distribution results in the abnormal transport are shown in Fig.
2C&F. The comparison between normal and abnormal transport shows that the
velocity magnitude decreases in the red dashed circle region due to the
reduced number of MTs, and this further leads to accumulation of the material
in this area.

Table 1: Simulation parameters utilized in computations Parameter | Description | Default value
---|---|---
$D_{\pm}$ | Diffusion coefficient of incoming and outgoing materials | 0.1
$k_{\pm}$ | Attachment rate to the MTs that transport materials in the positive (+) and negative (-) directions | 1.0
$k^{\prime}_{\pm}$ | Detachment rate from MTs for materials that move in the positive (+) and negative (-) directions | 0.1
$l_{\pm}$ | Density of MTs used for motor-assisted transport | 1.0
$\mu$ | Viscosity of the traffic flow | 0.1
$\lambda_{i}$ | Degree of loading at inlet end | 2.0
$\lambda_{o}$ | Degree of loading at outlet end | 2.0
$n_{i}$ | Boundary value of $n_{0}$ at inlet end | 1.0
$n_{o}$ | Boundary value of $n_{0}$ at outlet end | 0.0
$\alpha$ | Penalty parameter for the cost to control high concentration gradient | 1.0
$\beta$ | Penalty parameter for the cost of control force $f_{\pm}$ | 1.0

As shown in Fig. 2G-I and Fig. S1, we also perform parameter analysis using
the single pipe geometry to study the influence of simulation parameters on
the material distribution results. In particular, we focus on three parameters
that may have significant effect when dealing with traffic jam caused by the
reduction of MTs. The values selected for these parameters are displayed in
Table 2. We assume the active regulation from neuron is less dominant than
natural transport via diffusion or MTs, and thus select two smaller values for
$\alpha$ and $\beta$ compared to the default values in Table 1. Regarding the
value selection of $k/k^{\prime}$, we refer to the values utilized in [71] and
ensure the selected values stay within a biologically realistic range. Fig. 2G
shows the effect of the penalty parameter of the concentration gradient cost,
$\alpha$, on the concentration distribution. One can see that the decrease of
$\alpha$ leads to a severer material accumulation around the region with
reduced MTs in the single pipe geometry. We also find that the concentration
gradient becomes larger around the traffic jam region, which indicates that
there is less control over the concentration gradient due to the decrease of
$\alpha$.

Table 2: Value selection for parameter study Parameter | Value selection
---|---
$\alpha$ | 1, 0.1, 0.01
$\beta$ | 1, 0.1, 0.01
$k_{\pm}/k^{\prime}_{\pm}$ | Fix $k_{\pm}=1.0\;s^{-1}$, let $k_{\pm}/k^{\prime}_{\pm}$ = 1, 10, 100

Fig. 2H is similar to Fig. 2G but shows the effect of the penalty parameter of
the control force, $\beta$, on the concentration distribution. We find similar
phenomena that the traffic jam gets worse when $\beta$ decreases. By comparing
Fig. 2G with 2H and Fig. S1A with S1B, we find $\beta$ has a greater influence
on the concentration than $\alpha$ when decreasing both parameters by the same
amount. Since $\beta$ affects the control force in Eq. 1d while $\alpha$
affects the concentration in Eq. 1b&1c, the result indicates that the
regulation of transport velocity on MTs is vital to achieve the optimal
material transport process in neuron.

Fig. 2I shows the effect of the ratio between the attachment rate and
detachment rate, $k/k^{\prime}$, on the material concentration. We find that
when $k/k^{\prime}$ increases, the location of maximum concentration moves
toward right, which indicates the decrease of detachment rate $k^{\prime}$
causes more material get attached to MTs and transport faster as expected in
[18]. However, the reduction of MTs slows down the motor-assisted transport on
MTs and results in worse traffic jam. Interestingly, when $k/k^{\prime}$
decreases from 10 to 1, we also observe a similar traffic jam phenomena. The
possible reason is that the increase of $k^{\prime}$ causes more material
transported via free diffusion. Although free diffusion helps to transport the
material farther along the branch, the slow diffusion speed limits its ability
to mitigate the traffic jam caused by the reduction of MTs.

Figure 2: Simulation of material transport and parameter analysis in a single
pipe geometry. (A) The predefined velocity field $V_{+}$. Black arrow points
to the inlet of the pipe. The computed velocity field in (B) a healthy neuron
and (C) an abnormal neuron with reduced MTs in the red dashed circle region.
(D) Distribution of $l_{+}$ to model the traffic jam caused by the reduction
of MTs. Distribution of concentration in (E) a healthy neuron and (F) an
abnormal neuron with reduced MTs in the red dashed circle region. (G-I) The
concentration curve on the centerline of the single pipe affected by different
settings of (G) $\alpha$; (H) $\beta$; and (I) $k/k^{\prime}$. Unit for color
bars: (A-C) $\mu m/s$ and (E, F) $mol/\mu m^{3}$.

To account for morphological effect on the transport process, we simulate the
normal material transport and the abnormal transport with traffic jam in two
neuron tree structures as shown in Figs. 3&4. The predefined velocity fields
for both geometries are shown in Figs. 3A&4A. To quantitatively study the
influence of traffic jam on the material concentration among tree structures,
we also plot the concentration distribution curves along the centerline from
the inlet to each outlet of these two neurons. In each curve plot, we compare
the distribution between the normal transport and the abnormal transport with
traffic jam, as shown in Figs. 3E&4E. For both cases, we model traffic jam by
reducing the number of MTs ($l_{+}$) used for transport in the red dashed
circle regions. As a result, a sudden decrease of velocity (Figs. 3C&4C) and
material accumulation (Figs. 3E&4E) can be observed in these regions. By
observing the distribution curve of the outlets downstream the traffic jam
region (curve plots 1-4 of Fig. 3E and 3-8 of Fig. 4E we find that the reduced
number of MTs not only causes high concentration in the local region, but also
decreases the material concentration along the downstream of traffic jam
region. The distribution curves of the other outlets (curve plots 5 of Fig. 3E
and 1, 2, 9, 10 of Fig. 4E) demonstrate that more materials are transported to
these outlets to minimize the hazard of traffic jam. The result also shows
that materials rely on motor-assisted transport in longer branches of neurons
and the directional transport on MTs contributes significantly to the entire
transport process.

As shown in Fig. 3F-H and Fig. S2, we also perform parameter analysis on the
concentration distribution in the neuron tree structure. Similar to the
parameter analysis in single pipe geometry, we study the influence of three
parameters on the concentration distribution and the selected values are
listed in Table. 2. To quantitatively study the influence, we also plot and
compare the concentration curves on the centerline from inlet to outlet 2 of
the neuron tree. We obtain similar results as in single pipe geometry that the
decrease of $\alpha$ or $\beta$ leads to a severer material accumulation
around the region with reduced MTs, and $\beta$ shows greater effect than
$\alpha$ on the concentration distribution. In addition, we observe in Figs.
3E, 3F and S2 that when $\alpha$ or $\beta$ increases, more material is
transported to the bottom long branch to mitigate the traffic jam in other
branches. In Fig. 3H, we also find that the maximum concentration location
moves downstream slightly when $k/k^{\prime}$ increases, and either increasing
or decreasing $k/k^{\prime}$ intensifies the traffic jam.

Figure 3: Simulation of material transport and parameter analysis in a neuron
tree extracted from NMO_54504. (A) The predefined velocity field $V_{+}$.
Black arrow points to the inlet of the neuron tree. The computed velocity
field in (B) a healthy neuron and (C) an abnormal neuron with reduced MTs in
the red dashed circle region. Distribution of concentration in (D) a healthy
neuron and (E) an abnormal neuron. We also compare the concentration curve on
the centerline from the inlet to every outlet between normal and abnormal
transport in (E). The red dashed curve shows the centerline from the inlet to
one of the outlets and each outlet is indexed by a unique number. (F-H) The
concentration curve on the centerline from inlet to outlet 2 affected by
different settings of (F) $\alpha$; (G) $\beta$; and (H) $k/k^{\prime}$. Unit
for color bars: (A-C) $\mu m/s$ and (D, E) $mol/\mu m^{3}$. Figure 4:
Simulation of material transport in a neuron tree extracted from NMO_54499.
(A) The predefined velocity field $V_{+}$. Black arrow points to the inlet of
the material. The computed velocity field in (B) a healthy neuron and (C) an
abnormal neuron with reduced MTs in the red dashed circle region. Distribution
of concentration and the concentration curve on the centerline of the circled
region in (D) a healthy neuron and (E) an abnormal neuron. We also compare the
concentration curve on the centerline from the inlet to every outlet between
normal and abnormal transport in (E). The red dashed curve shows the
centerline from the inlet to one of the outlets and each outlet is indexed by
a unique number. Unit for color bars: (A-C) $\mu m/s$ and (D, E) $mol/\mu
m^{3}$.

### Simulation of traffic jam with MT swirls and local swelling

Recent studies have shown that the formation of MT swirls can lead to
accumulation of transported material and cause local swelling of neuron
geometries [72]. In our model, we modify the spatial distribution of $l_{\pm}$
and enlarge the radius of neuron in a local region to simulate the effect of
MT swirls and local swelling on the transport process. We explain the
simulation setting by using a straight pipe geometry with MT swirls and
swelling in the middle $L_{2}$ region, as shown in Fig. 5. We assume the
normal transport is unidirectional from left to right ($+$ direction, red
arrow in Fig. 5A). Due to the MT swirls in the middle region, the transport
path is extended by two segments: one segment reverses to transport the
material from right to left ($-$ direction, blue arrow in Fig. 5A) and the
other segment transports in the normal direction from left to right.
Therefore, we increase the values of $l_{+}$ and $l_{-}$ along the
longitudinal direction in the swelling region to describe the transport path
change caused by swirling. We also assume that the swirl direction is counter-
clockwise and assign an asymmetric distribution of $l_{\pm}$ on the cross-
section in the swelling region. In particular, $l_{+}$ is higher on the bottom
of cross-section while $l_{-}$ is higher on the top of cross-section, as shown
in Fig. 5A. We perform simulation with the new parameter setting and compare
with the results of normal transport in the same geometry. The velocity field
and concentration distribution of normal and abnormal transport are compared
in Fig. 5B&C, respectively. The decrease of velocity and material accumulation
can be observed in the swollen region. We also find that the velocity
magnitude is not symmetric anymore due to the MT swirls in abnormal transport.
In Fig. 5D, we plot the velocity streamline with concentration distribution in
the zoomed-in swollen region for both normal and abnormal transport. Compared
to the uniform velocity streamline in normal transport, the velocity displays
vortex pattern in the abnormal transport, which reflects a longer transport
distance due to MT swirls. We also find that the swirl of velocity streamline
usually happens in the high concentration region, which implies that the
material accumulation is caused by the vortex-shape velocity field.

As shown in Fig. 6, we then apply the same approach to simulate the normal and
abnormal transport with MT swirls in two neuron tree structures with local
swelling. The swelling is introduced by increasing the skeleton radius in the
red dashed circle regions. We also assume a counter-clockwise MT swirl in
these swollen regions and modify the distribution of $l_{\pm}$ accordingly.
For each model, we simulate the abnormal transport process due to MT swirls to
obtain velocity field and concentration distribution results and compare with
the result of normal transport in the same geometry. By comparing Fig. 6A&C
with Fig. 6B&D, we find the velocity magnitude decreases and material
accumulates in the swollen region. In other branches that are not downstream
the swollen region, the material concentration also increases to mitigate the
traffic jam in the swollen region. In addition, similar to the results in
straight pipe with swelling geometry (Fig. 5D), we also observe that the
velocity streamline with vortex pattern matches with the high concentration
region (Fig. 6B&D). These results illustrate that the MT swirls lead to the
circular transport velocity field in a local region which not only extends the
transport distance but also traps the material and causes traffic jam.

Figure 5: Simulation of material transport in a straight pipe with swelling in
the middle region. (A) The simulation setting for modeling MT swirls. The red
and blue arrows show the transport path along swirly MTs. Due to the MT swirls
in the $L_{2}$ region, both $l_{+}$ and $l_{-}$ are increased along centerline
and their distributions on cross-section are also modified. (B, C) The
computed velocity field and concentration distribution in the swollen
geometry. (D) The velocity streamline and concentration distribution in the
swollen region. Different color maps are used to distinguish between velocity
and concentration. Unit for color bars: Concentration: $mol/\mu m^{3}$;
Velocity: $\mu m/s$. Figure 6: Simulation of material transport in neuron
trees extracted from (A, B) NMO_54504 and (C, D) NMO_54499 with swelling in
the red circle regions. The first column shows the computed velocity field and
black arrow points to the inlet of the material. The second column shows the
concentration distribution. The last column shows the velocity streamline and
concentration distribution in the swollen region. Different color maps are
used to distinguish between velocity and concentration. Unit for color bars:
Concentration: $mol/\mu m^{3}$; Velocity: $\mu m/s$.

## Discussion

In this paper, we develop a PDE-constrained optimization model to simulate
material transport control in neurons. Using our simulation, we examine both
normal and abnormal transport processes in different geometries and discover
several spatial patterns of the transport process. Our results show the
formation of traffic jams due to the reduction of MTs and MT swirls in the
local region. We also observe how the traffic jam affects the spatial patterns
of transport velocities that in turn drives the transported materials
distributed distinctly in different regions of neurite networks to mitigate
traffic jam. By solving the proposed new optimization problem, we build a more
realistic transport model for neurons by including active traffic regulation.
The model is successfully applied to complex 2D neuron geometries and provides
key insights into how neuron mediates the material transport inside its
complex geometry.

Our study shows that MTs have a major impact on the material transport
velocity and further affect the material concentration distribution. As shown
in Fig. 2, the reduction of MTs in the middle of the single pipe slows down
the transport velocity downstream and leads to traffic jam in the middle
region. When the neuron has more branches in its geometry (Figs. 3&4), the
reduction of MTs in one branch has a similar influence on the transport
downstream the branch. However, we observe an increase in transport velocity
and material concentration in other branches, indicating that the active
regulation from neuron takes effect to avoid traffic jams. In addition, we
perform parameter analysis to study the influence of different simulation
parameters on the material concentration distribution. The ratio between the
attachment rate $k$ and detachment rate $k^{\prime}$ affects the amount of
material transported via MTs or free diffusion. This will affect the overall
transport speed and material distribution due to the different transport
behaviour between motor-assisted transport and free diffusion. The penalty
parameters $\alpha$ and $\beta$ affect the ability of neuron to handle traffic
jams. $\beta$ has a greater influence on the traffic regulation compared to
$\alpha$ since it directly affects the transport velocity on MTs, this again
verifies the vital role of MTs during the intracellular transport process. Our
model can also model the influence of diverse neuron topologies on material
distribution. For the transport in healthy neurons (Figs. 3-4B&D), the
magnitude of transport velocity is different among branches due to the
asymmetric geometry. The different velocity magnitude further contributes to
the distinct material concentration in different branches. In particular, we
find that shorter branches tend to have faster transport speed and higher
material concentration, which may result from the high demand of materials for
their growth.

Our study also successfully simulates and provides reasonable explanation on
the traffic jam caused by MT swirls. We assume the counter-clockwise MT swirls
exist in a local region of neuron geometry which cause traffic jam and
geometry swelling. The spatial distribution of MT density ($l_{\pm}$) and
neuron geometry are modified accordingly to model this phenomena. We compare
the simulation result of abnormal transport on swirly MTs with normal
transport and find that MT swirls have severe impact on the transport velocity
field. Compared to the uniform velocity streamline in normal transport, the
abnormal transport exhibits a streamline with counter-clockwise vortex pattern
(Figs. 5D, 6B&D), which is caused by the counter-clockwise MT swirls. This
circular streamline not only extends the transport distance but also traps the
material in the local region, and therefore explains why high concentration
region matches with the circular streamline pattern.

Our study develops an IGA solver (available at
https://github.com/truthlive/NeuronTransportOptimization) for solving the PDE-
CO problem in complex neuron geometries. Specifically, we adopt the skeleton-
based sweeping method [32, 38] for mesh generation to represent the tree
structures of neuron geometry. Given the geometry information of neurons, our
method automatically reconstructs 2D network geometry with high accuracy and
high order of continuity for IGA computation. Our automatic IGA optimization
solver provides an efficient computation tool for studies of material
transport regulation in complex neurite networks. The current 2D solver can be
easily generalized to 3D and it is also extensible to solve other PDE-CO
models of cellular processes in complex neurite network geometry.

Our study has its limitations, which we are addressing in the ongoing work. In
the current model, we only consider the influence of traffic jams on the
material concentration but neglect its effect on the deformation of neuron
geometries. In addition, although IGA offers great advantages in accurately
simulating material transport control in complex neuron geometries, the
computational cost of simulating transport in large-scale neurite networks
remains very expensive, which limits its biomedical application. To improve
the computational efficiency of our model, we will adopt deep learning
techniques to build fast and accurate surrogate models [73, 50]. Despite these
limitations, our simulation directly shows how the traffic jam is formed in
neurons and how neurons could control material traffic to avoid traffic jams.
The simulation results provide references to further answer the question of
how neurons deliver the right material to the right destination in a balanced
manner in their complex neurite networks and how the transport may be affected
by disease conditions.

## Code and data availability

The source code for our model and all input data are available for download
from a public software repository located at
https://github.com/truthlive/NeuronTransportOptimization. All data generated
during this study can be reconstructed by running the source code.

## Acknowledgememts

The authors acknowledge the support of NSF grants CMMI-1953323 and
CBET-1804929. This work used the Extreme Science and Engineering Discovery
Environment (XSEDE), which is supported by National Science Foundation grant
number ACI-1548562. Specifically, it used the Bridges-2 system, which is
supported by NSF award number ACI-1928147, at the Pittsburgh Supercomputing
Center (PSC).

## References

  * [1] Idan Segev and Michael London.  Untangling dendrites with quantitative models.  Science, 290(5492):744–750, 2000.
  * [2] Sharon A Swanger and Gary J Bassell.  Dendritic protein synthesis in the normal and diseased brain.  Neuroscience, 232:106–127, 2013.
  * [3] Kurt J De Vos, Andrew J Grierson, Steven Ackerley, and Christopher CJ Miller.  Role of axonal transport in neurodegenerative diseases.  Nature Review Neuroscience, 31:151–173, 2008.
  * [4] Shermali Gunawardena and Lawrence SB Goldstein.  Polyglutamine diseases and transport problems: deadly traffic jams on neuronal highways.  Archives of Neurology, 62(1):46–51, 2005.
  * [5] Stéphanie Millecamps and Jean-Pierre Julien.  Axonal transport deficits and neurodegenerative diseases.  Nature Reviews Neuroscience, 14(3):161–176, 2013.
  * [6] Natalia L Kononenko, Gala A Claßen, Marijn Kuijpers, Dmytro Puchkov, Tanja Maritzen, Aleksandra Tempes, Anna R Malik, Agnieszka Skalecka, Sujoy Bera, Jacek Jaworski, et al.  Retrograde transport of TrkB-containing autophagosomes via the adaptor AP-2 mediates neuronal complexity and prevents neurodegeneration.  Nature Communications, 8(1):1–16, 2017.
  * [7] Hongyu Zhang, Chunlei Zhang, Jean Vincent, Diana Zala, Caroline Benstaali, Matthieu Sainlos, Dolors Grillo-Bosch, Sophie Daburon, Françoise Coussen, Yoon Cho, et al.  Modulation of AMPA receptor surface diffusion restores hippocampal plasticity and memory in Huntington’s disease models.  Nature Communications, 9(1):1–16, 2018.
  * [8] Ronald D Vale.  The molecular motor toolbox for intracellular transport.  Cell, 112(4):467–480, 2003.
  * [9] Nobutaka Hirokawa, Shinsuke Niwa, and Yosuke Tanaka.  Molecular motors in neurons: transport mechanisms and roles in brain function, development, and disease.  Neuron, 68(4):610–638, nov 2010.
  * [10] Mariella AM Franker and Casper C Hoogenraad.  Microtubule-based transport–basic mechanisms, traffic rules and role in neurological pathogenesis.  Journal of Cell Science, 126(11):2319–2329, 2013.
  * [11] MGL Van den Heuvel, MP De Graaff, and C Dekker.  Microtubule curvatures under perpendicular electric forces reveal a low persistence length.  Proceedings of the National Academy of Sciences, 105(23):7941–7946, 2008.
  * [12] Shoichiro Tsukita and Harunori Ishikawa.  The cytoskeleton in myelinated axons: serial section study.  Biomedical Research, 2(4):424–437, 1981.
  * [13] Bl J Schnapp and TS Reese.  Cytoplasmic structure in rapid-frozen axons.  The Journal of Cell Biology, 94(3):667–669, 1982.
  * [14] Nobutaka Hirokawa and Yasuko Noda.  Intracellular transport and kinesin superfamily proteins, KIFs: structure, function, and dynamics.  Physiological Reviews, 88(3):1089–1118, 2008.
  * [15] Scott R May, Amir M Ashique, Mattias Karlen, Baolin Wang, Yiguo Shen, Kostantinos Zarbalis, Jeremy Reiter, Johan Ericson, and Andrew S Peterson.  Loss of the retrograde motor for IFT disrupts localization of Smo to cilia and prevents the expression of both activator and repressor functions of Gli.  Developmental Biology, 287(2):378–389, 2005.
  * [16] Robert Day Allen, Janis Metuzals, Ichiji Tasaki, Scott T Brady, and Susan P Gilbert.  Fast axonal transport in squid giant axon.  Science, 218(4577):1127–1129, 1982.
  * [17] Michael C Reed, Stephanos Venakides, and Jacob Joseph Blum.  Approximate traveling waves in linear reaction-hyperbolic equations.  SIAM Journal on Applied Mathematics, 50(1):167–180, 1990.
  * [18] DA Smith and RM Simmons.  Models of motor-assisted transport of intracellular particles.  Biophysical Journal, 80(1):45–68, 2001.
  * [19] Avner Friedman and Gheorghe Craciun.  A model of intracellular transport of particles in an axon.  Journal of Mathematical Biology, 51(2):217–246, 2005.
  * [20] Elizabeth A. Brooks.  Probabilistic methods for a linear reaction-hyperbolic system with constant coefficients.  The Annals of Applied Probability, 9(3):719 – 731, 1999.
  * [21] Lea Popovic, Scott A McKinley, and Michael C Reed.  A stochastic compartmental model for fast axonal transport.  SIAM Journal on Applied Mathematics, 71(4):1531–1556, 2011.
  * [22] Jay M Newby and Paul C Bressloff.  Quasi-steady state reduction of molecular motor-based models of directed intermittent search.  Bulletin of Mathematical Biology, 72(7):1840–1866, 2010.
  * [23] Paul C Bressloff and Jay M Newby.  Stochastic hybrid model of spontaneous dendritic NMDA spikes.  Physical Biology, 11(1):016006, 2014.
  * [24] Chuan Xue, Blerta Shtylla, and Anthony Brown.  A stochastic multiscale model that explains the segregation of axonal microtubules and neurofilaments in neurological diseases.  PLoS Comput Biol, 11(8):e1004406, 2015.
  * [25] Michiel Bertsch, Bruno Franchi, Norina Marcello, Maria Carla Tesi, and Andrea Tosin.  Alzheimer’s disease: a mathematical model for onset and progression.  Mathematical Medicine and Biology: A Journal of the IMA, 34(2):193–214, 2017.
  * [26] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  Computer Methods in Applied Mechanics and Engineering, 194(39):4135 – 4195, 2005.
  * [27] Les Piegl and Wayne Tiller.  The NURBS Book.  Springer Science & Business Media, 2012.
  * [28] DJ Benson, Y Bazilevs, Ming-Chen Hsu, and Thomas JR Hughes.  Isogeometric shell analysis: the Reissner–Mindlin shell.  Computer Methods in Applied Mechanics and Engineering, 199(5-8):276–289, 2010.
  * [29] Hugo Casquero, Lei Liu, Yongjie Zhang, Alessandro Reali, Josef Kiendl, and Hector Gomez.  Arbitrary-degree T-splines for isogeometric analysis of fully nonlinear Kirchhoff–Love shells.  Computer-Aided Design, 82:140–153, 2017.
  * [30] Hugo Casquero, Xiaodong Wei, Deepesh Toshniwal, Angran Li, Thomas JR Hughes, Josef Kiendl, and Yongjie Jessica Zhang.  Seamless integration of design and Kirchhoff–Love shell analysis using analysis-suitable unstructured T-splines.  Computer Methods in Applied Mechanics and Engineering, 360:112765, 2020.
  * [31] Xiaodong Wei, Xin Li, Kuanren Qian, Thomas J. R. Hughes, Yongjie Jessica Zhang, and Hugo Casquero.  Analysis-suitable unstructured T-splines: multiple extraordinary points per face, 2021.
  * [32] Yongjie Zhang, Yuri Bazilevs, Samrat Goswami, Chandrajit L Bajaj, and Thomas JR Hughes.  Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow.  Computer Methods in Applied Mechanics and Engineering, 196(29-30):2943–2959, 2007.
  * [33] Yongjie Zhang, Xinghua Liang, Jun Ma, Yiming Jing, Matthew J Gonzales, Christopher Villongco, Adarsh Krishnamurthy, Lawrence R Frank, Vishal Nigam, Paul Stark, et al.  An atlas-based geometry pipeline for cardiac Hermite model construction and diffusion tensor reorientation.  Medical Image Analysis, 16(6):1130–1141, 2012.
  * [34] Yongjie Zhang.  Challenges and advances in image-based geometric modeling and mesh generation.  In Yongjie Jessica Zhang, editor, Image-Based Geometric Modeling and Mesh Generation, pages 1–10. Springer, 2013.
  * [35] Benjamin Urick, Travis M Sanders, Shaolie S Hossain, Yongjie J Zhang, and Thomas JR Hughes.  Review of patient-specific vascular modeling: template-based isogeometric framework and the case for CAD.  Archives of Computational Methods in Engineering, 26(2):381–404, 2019.
  * [36] Yuxuan Yu, Yongjie Jessica Zhang, Kenji Takizawa, Tayfun E Tezduyar, and Takafumi Sasaki.  Anatomically realistic lumen motion representation in patient-specific space–time isogeometric flow analysis of coronary arteries with time-dependent medical-image data.  Computational Mechanics, 65(2):395–404, 2020.
  * [37] Yongjie Jessica Zhang.  Geometric Modeling and Mesh Generation from Scanned Images, volume 6.  CRC Press, 2016.
  * [38] Angran Li, Xiaoqi Chai, Ge Yang, and Yongjie Jessica Zhang.  An isogeometric analysis computational platform for material transport simulation in complex neurite networks.  Molecular & Cellular Biomechanics, 16(2):123–140, 2019.
  * [39] Aishwarya Pawar and Yongjie Jessica Zhang.  NeuronSeg_BACH: automated neuron segmentation using B-Spline based active contour and hyperelastic regularization.  Communications in Computational Physics, 28(3):1219–1244, 2020\.
  * [40] Yuri Bazilevs, Victor M Calo, Yongjie Zhang, and Thomas JR Hughes.  Isogeometric fluid–structure interaction analysis with applications to arterial blood flow.  Computational Mechanics, 38(4-5):310–322, 2006.
  * [41] Hugo Casquero, Lei Liu, Carles Bona-Casas, Yongjie Zhang, and Hector Gomez.  A hybrid variational-collocation immersed method for fluid-structure interaction using unstructured T-splines.  International Journal for Numerical Methods in Engineering, 105(11):855–880, 2016.
  * [42] Hugo Casquero, Yongjie Jessica Zhang, Carles Bona-Casas, Lisandro Dalcin, and Hector Gomez.  Non-body-fitted fluid–structure interaction: divergence-conforming B-splines, fully-implicit dynamics, and variational formulation.  Journal of Computational Physics, 374:625–653, 2018.
  * [43] Hugo Casquero, Carles Bona-Casas, Deepesh Toshniwal, Thomas JR Hughes, Hector Gomez, and Yongjie Jessica Zhang.  The divergence-conforming immersed boundary method: application to vesicle and capsule dynamics.  Journal of Computational Physics, 425:109872, 2021.
  * [44] Yuxuan Yu, Xiaodong Wei, Angran Li, Jialei Ginny Liu, Jeffrey He, and Yongjie Jessica Zhang.  HexGen and Hex2Spline: polycube-based hexahedral mesh generation and spline modeling for isogeometric analysis applications in LS-DYNA.  In Springer INdAM Serie: Proceedings of INdAM Workshop "Geometric Challenges in Isogeometric Analysis". Springer, 2020.
  * [45] Yuxuan Yu, Jialei Ginny Liu, and Yongjie Jessica Zhang.  HexDom: polycube-based hexahedral-dominant mesh generation.  In Kenneth Morgan Rebén Sevilla, Simona Perotto, editor, The Edited Volume of Mesh Generation and Adaptation: Cutting-Edge Techniques for the 60th Birthday of Oubay Hassan. SEMA-SIMAI Springer Series. Springer, 2021\.
  * [46] Xiaodong Wei, Yongjie Zhang, Lei Liu, and Thomas JR Hughes.  Truncated T-splines: fundamentals and methods.  Computer Methods in Applied Mechanics and Engineering, 316:349–372, 2017.
  * [47] Xiaodong Wei, Yongjie Jessica Zhang, and Thomas JR Hughes.  Truncated hierarchical tricubic $C^{0}$ spline construction on unstructured hexahedral meshes for isogeometric analysis applications.  Computers & Mathematics with Applications, 74(9):2203–2220, 2017\.
  * [48] Xiaodong Wei, Yongjie Jessica Zhang, Deepesh Toshniwal, Hendrik Speleers, Xin Li, Carla Manni, John A Evans, and Thomas JR Hughes.  Blended B-spline construction on unstructured quadrilateral and hexahedral meshes with optimal convergence rates in isogeometric analysis.  Computer Methods in Applied Mechanics and Engineering, 341:609–639, 2018.
  * [49] Xin Li, Xiaodong Wei, and Yongjie Jessica Zhang.  Hybrid non-uniform recursive subdivision with improved convergence rates.  Computer Methods in Applied Mechanics and Engineering, 352:606–624, 2019.
  * [50] Angran Li, Amir Barati Farimani, and Yongjie Jessica Zhang.  Deep learning of material transport in complex neurite networks.  Scientific Reports, 11(1):1–13, 2021.
  * [51] Jeff Borggaard and John Burns.  A PDE sensitivity equation method for optimal aerodynamic design.  Journal of Computational Physics, 136(2):366–384, 1997.
  * [52] Christian Brandenburg, Florian Lindemann, Michael Ulbrich, and Stefan Ulbrich.  Advanced numerical methods for PDE constrained optimization with application to optimal design in Navier Stokes flow.  In Constrained optimization and optimal control for partial differential equations, pages 257–275. Springer, 2012.
  * [53] Michael Hinze, René Pinnau, Michael Ulbrich, and Stefan Ulbrich.  Optimization with PDE Constraints, volume 23.  Springer Science & Business Media, 2008.
  * [54] Lorenz T Biegler, Omar Ghattas, Matthias Heinkenschloss, David Keyes, and Bart van Bloemen Waanders.  Real-time PDE-constrained Optimization.  SIAM, 2007.
  * [55] Roland Herzog and Karl Kunisch.  Algorithms for PDE-constrained optimization.  GAMM-Mitteilungen, 33(2):163–176, 2010.
  * [56] Tyrone Rees, H Sue Dollar, and Andrew J Wathen.  Optimal solvers for PDE-constrained optimization.  SIAM Journal on Scientific Computing, 32(1):271–298, 2010.
  * [57] Eldad Haber, Matthias Chung, and Felix Herrmann.  An effective method for parameter estimation with PDE constraints with multiple right-hand sides.  SIAM Journal on Optimization, 22(3):739–757, 2012.
  * [58] Xiaolei Xun, Jiguo Cao, Bani Mallick, Arnab Maity, and Raymond J Carroll.  Parameter estimation of partial differential equation models.  Journal of the American Statistical Association, 108(503):1009–1020, 2013.
  * [59] Cosmina Hogea, Christos Davatzikos, and George Biros.  An image-driven parameter estimation problem for a reaction–diffusion glioma growth model with mass effects.  Journal of Mathematical Biology, 56(6):793–825, 2008.
  * [60] Hyun Keol Kim, Molly Flexman, Darrell J Yamashiro, Jessica J Kandel, and Andreas H Hielscher.  Transport-theory based multispectral imaging with PDE-constrained optimization.  In Optical Tomography and Spectroscopy of Tissue IX, volume 7896, page 78960W. International Society for Optics and Photonics, 2011.
  * [61] Alessandro Melani.  Adjoint-based parameter estimation in human vascular one dimensional models.  PhD thesis, Mathematical Models and Methods in Engineering, Department of Mathematics, Politecnico di Milano, 2013.
  * [62] Damián Alejandro Knopoff, Damián R Fernández, Germán Ariel Torres, and Cristina Vilma Turner.  Adjoint method for a tumor growth PDE-constrained optimization problem.  Computers & Mathematics with Applications, 66(6):1104–1119, 2013\.
  * [63] Andrés Agustín Ignacio Quiroga, Damián Fernández, Germán Ariel Torres, and Cristina Vilma Turner.  Adjoint method for a tumor invasion PDE-constrained optimization problem in 2D using adaptive finite element method.  Applied Mathematics and Computation, 270:358–368, 2015.
  * [64] Martin Stoll and Andy Wathen.  All-at-once solution of time-dependent Stokes control.  Journal of Computational Physics, 232(1):498–515, 2013.
  * [65] Fikriye Yılmaz and Bülent Karasözen.  An all-at-once approach for the optimal control of the unsteady Burgers equation.  Journal of Computational and Applied Mathematics, 259:771–779, 2014\.
  * [66] Youcef Saad and Martin H Schultz.  GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986.
  * [67] Shrirang Abhyankar, Jed Brown, Emil M. Constantinescu, Debojyoti Ghosh, Barry F. Smith, and Hong Zhang.  PETSc/TS: a modern scalable ode/dae solver library, 2018.
  * [68] Giorgio A Ascoli, Duncan E Donohue, and Maryam Halavi.  NeuroMorpho.Org: a central resource for neuronal morphologies.  Journal of Neuroscience, 27(35):9247–9251, 2007.
  * [69] John Towns, Timothy Cockerill, Maytal Dahan, Ian Foster, Kelly Gaither, Andrew Grimshaw, Victor Hazlewood, Scott Lathrop, Dave Lifka, Gregory D Peterson, et al.  XSEDE: accelerating scientific discovery.  Computing in Science & Engineering, 16(5):62–74, 2014.
  * [70] Nicholas A Nystrom, Michael J Levine, Ralph Z Roskies, and J Ray Scott.  Bridges: a uniquely flexible HPC resource for new communities and data analytics.  In Proceedings of the 2015 XSEDE Conference: Scientific Advancements Enabled by Enhanced Cyberinfrastructure, pages 1–8, 2015.
  * [71] AV Kuznetsov and AA Avramenko.  A macroscopic model of traffic jams in axons.  Mathematical Biosciences, 218(2):142–152, 2009.
  * [72] Or A Shemesh, Hadas Erez, Irith Ginzburg, and Micha E Spira.  Tau-induced traffic jams reflect organelles accumulation at points of microtubule polar mismatching.  Traffic, 9(4):458–471, 2008.
  * [73] Angran Li, Ruijia Chen, Amir Barati Farimani, and Yongjie Jessica Zhang.  Reaction diffusion system prediction based on convolutional neural network.  Scientific Reports, 10(1):1–9, 2020.

## Supplementary Information

Figure S1: Parameter analysis in the single pipe geometry. The distribution of
concentration affected by different settings of (A) penalty parameter $\alpha$
for the cost to control high concentration gradient; (B) penalty parameter
$\beta$ for the cost of control force; and (C) the ratio between attachment
rate $k$ and detachment rate $k^{\prime}$. Unit for color bars: $mol/\mu
m^{3}$. Figure S2: Parameter analysis in the neuron tree extracted from
NMO_54504. The distribution of concentration affected by different settings of
(A) penalty parameter $\alpha$ for the cost to control high concentration
gradient; (B) penalty parameter $\beta$ for the cost of control force; and (C)
the ratio between attachment rate $k$ and detachment rate $k^{\prime}$. Unit
for color bars: $mol/\mu m^{3}$.