Threshold dynamics of a stochastic SIRS epidemic model with transfer from infected individuals to susceptible individuals and log-normal Ornstein-Uhlenbeck process

1.   Introduction

This paper is an invitation to consider a class of what I may call experimental Lie algebras and special functions associated with them. The hypothesis is that such explorations may be useful in understanding the phenomena from which they arise naturally, namely the formation of biological organs. For illustrative purposes, we shall focus on biological growth patterns as viewed at the mesoscale: the level where continuum physics remains useful, with spatial resolution well above the scale of individual cells. For all our examples and data we shall specialize to the growth of the brain and of the cortical convolutions that have fascinated so many people. The approach however is much broader than that application as should become clear later. Physical models are here understood to address mechanical forces, leaving the biochemical and genetic substrates of growth to act as external forces. Such models of biological growth have by now acquired a tradition, starting perhaps with D'Arcy Wentworth Thompson ^[1].

In this Introduction, the outline of the rest of this paper is described, and remarks on the literature of the subjects considered are mentioned. The theme of the special issue for which this paper is submitted is the mathematics of symmetry groups and equivalence due in main to Sophus Lie and Elie Cartan as applied to equations (primarily differential and integro-differential) describing biological phenomena. We therefore begin the next section (section 2) with brief references to the mathematical modeling of brain and cortical growth, and explore the very simplest equations that model growth. We describe the generalization of the method of symmetries required to give the generators of growth, and the Lie algebra thereby generated. However, we find even in this simplest model that formulas which would result seem too complex for implementation and use in the near term. Standard methods of Lie symmetries applied to equations may not be very useful in this area where the equations and constraints that define growth patterns are more complicated.

Section 3 offers another approach appropriate to this century of big data. Either with computer models of growth or with data collected from brain imaging, we can calculate numerical vector fields that generate growth. We provide examples both from our own work with brain images, and from other modern work that is available to the research community. These provide several different datasets from which generating vector fields of different dimensionalities may be obtained. Although the Lie algebra so generated can be of any, including formally infinite, dimensionality, the elements of the algebra are generated from a finite, even small, number of generators.

Section 4 is our preliminary exploration for how we might use such generators of Lie algebras. It is well known that Lie algebras (and groups) have special functions associated with them. The algebras we are concerned with are generalizations of nilpotent algebras such as the Heisenberg algebra. We describe several paths to obtaining special functions associated with the algebras in question, including infinite dimensional ones, and illustrate the lowest order special functions for particular generators, though it will be clear how the process can be carried out in general.

Models and data, Lie symmetries and special functions

In the rest of this introduction, we provide a guide to some literature on the topics we bring together in the rest of the paper, namely those mentioned in this subsection title, but primarily on the brain and models for its growth. The literature on computer models and related pertinent information on brain growth is truly vast, and a selection of these may be found here: ^{[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]}.

The normal gross morphology of the brain has been described in ^{[37,38,39,40,41,42,43]}. Older work on growth patterns ^[44,45] and neuronal migration: ^[46,47,48] continue to be developed. The morphology of the brain has always been a fascinating topic, and some of the more absurd conclusions drawn in relation to Einstein's brain may be found in ^[49,50]. Brain growth (ontogeny) is not the only area in which morphological patterns are interesting: they play a role in phylogeny (evolution) as well ^{[51,52,53,54]}. A discussion on the evolution of a particular sulcus called the lunate led to a controversy that may be traced in ^{[51,55,56,57,58,59,60,61,62]}. The gyrification or the folding of the brain is discussed in ^{[52,62,63,64]}. Cortical asymmetries ^[65,66,67] as well as sex differences ^[68] remain of interest as citation searches of the quoted references show. Finally, comparative morphology between species ^[69,70,71] deserves to be explored more fully. In fact, computer models of brain growth have mostly focused on the period from conception to gestation of the human, and more comprehensive models that integrate other species as well would be welcome.

The dream of a reductionist theory would be to construct the form of the brain from the gene expressions, and as can be surmised, much exciting work is being done on spatial patterns of gene expression, genetic control of cortical development and first-principles chemical kinetics ^{[72,73,74,75,76,77,78,79,80,81,82,83]}. In the context of cellular growth processes, morphogens and associated chemotaxis are extensively discussed in ^[84], which cites earlier work including ^[85,86,87]; see also ^[88,89] for an analysis of the chemical reactions and diffusions that generate and sustain morphogen gradients. However, "quite possibly there is no obligatory relationship between functional systems in the brain and global patterns of gene expression" ^[90], so that models restricted to the coarser scale of visible morphology remain of great interest. The radial patterns of migration of neurons from the ventricles have long been projected to be the cause of the tubular patterns that are so characteristic of the surface of the brain. A recent summary of such issues may be found in the thesis ^[91]. Much work still remains to synthesize the finer spatial scale at which gene expression occurs with the coarser scale of visible cortical folding patterns driven by mechanical constraints.

Patterns of folding and of sulci have analogs in entirely different fields, from the folding of dough ^[92,93,94] to geological folds ^[95,96,97]. Reaction-diffusion equations in two-dimensional space with certain coefficients have been demonstrated to give rise to labyrinthine tubular patterns ^[98], and this phenomenon has been invoked to explain certain features of cortical topography ^[99]. This planar theory may be relevant for three-dimensional tubular patterns because, generically, a section of an axisymmetric shell behaves similarly to any other section. We shall refer later to a more recent computer model using such a basis for describing the cortex. Finally, it is interesting to compare the pattern of sulci and their branching patterns such as found in ^[100] and to see if there are any "laws" analogous to those found in river meanders ^{[101,102,103]}. We now describe approaches to the three threads mentioned.

Mathematical models of growth were pioneered, as far as we are aware, by Richard Skalak and his collaborators and students (see ^[104] and references therein). They took as their starting point the continuum theory of dislocations. The chapter on dislocations and on disclinations in liquid crystals in the classical treatise on elasticity ^[105] can serve as a ready reference for the basic theory, while a full account of the continuum theory may be found in ^[106], where some of the history is given. In the last decade of the twentieth century, people used these concepts to describe growth in the continuum. Following the above reference, a few of a large number of papers utilizing this approach are ^{[107,108,109,110]}. Data on growth will be mentioned in section 3.

There are a large number of references on Lie symmetries. We have basically relied on the pedagogical masterpiece ^[111], in particular, several of the formulas in Chapter $5$ of that book. It should be understood that although we will not make repeated citations to ^[111], it, or other treatises on the Lie approach to differential equations, can be consulted for an explanation for all the material in section 2.2. Olver's later book ^[112] written at a higher mathematical level is not used here.

Finally, the papers on Lie algebras and special functions that we have consulted are mentioned in section 4. There is an encylopaedic three-volume series on the Representations of Lie Groups and Special Functions, followed by a fourth one on subsequent developments, by Vilenkin and Klimyk, but we have not consulted any of these.

2.   Lie algebras from growth-model equations

We first describe two models of cortical and brain growth, and then focus on the one that is more suited to illustrate how we might use the approach of Lie symmetries to attempt to discover the invariants and generators of growth. Section 2.2 then develops the equations of the symmetries. Some details on the formulas are deferred to Appendix A.

2.1.   Models

We shall refer only to two computer models of brain growth, in addition to our own older one. A three dimensional model that shows realistic patterns is described in ^[33] with the authors' previous relevant work referenced there. For the purposes of developing the first mathematical model, we will begin with a lower dimensional model. Below we reproduce pictorial representations of three mutually orthogonal sections of a human brain from our work creating an electronic brain atlas ^[113]. The sectioning planes are canonical, and are axial or horizontal for vertically standing human, sagittal, and coronal or frontal. These images are from the famous original paper atlases of Talairach and Tournoux, though here taken from the electronic version.

The cortical contours depicted in these sections are quite characteristic. Such cortical contours are the subject of a thesis ^[91] that develops a detailed model with the thickness of the cortex included in a planar model. The model is entirely numerical, and provides information useful for the purposes of this paper, as we shall describe in the next section. However, we turn to an older model that is less realistic but which expresses the key factors affecting cortical convolutions in a variational principle. This model was described in ^[114] and has the virtue for purposes of mathematical simplicity that it constructs an energy on a planar curve which is meant to be a cortical contour such as exhibited in the figure. The energy is minimized to obtain or predict the shape of the curve. The configuration of the curve is given by the position vector ${\bf{u}} (s) = \left(x(s), y(s)\right)$ of points on it. We will use the fact that the "coordinate" $s$ has exactly the same range before and after a transformation which alters the function ${\bf{u}}(s)$. Thus we can compare the values of the transformed position of the curve with the original at the same value of $s$. Further, as for the simpler case of geodesics for a given metric on a manifold, the flow is on the "jet bundle", that is, the curve is generated by minimizing an integral:

where $l(s)$ is the function that drives the growth. It is the rest length of the curve per unit $s$, along the curve. ${\bf{u}}_{J}$ denotes the $J^{\prime\rm{th}}$ derivative of ${\bf{u}}$: we omit the subscript $0$ when denoting ${\bf{u}}$ itself. The energy density $E$ consists of several terms describing physical effects: stretching, bending, skull roof and floor constraints, and those of white matter cabling and shear. The resulting energy density expressions involve differential functions of up to the second order so that stretching and bending are included. The behavior of white matter led us to introduce a term in $E$ that is an integral. In other words, $E$ is itself a functional. Thus as Figure 2(a) shows, the resulting thick curve is meant to be the cortical contour. The vault of the skull is represented by a semicircle, while the deep brain structures constraining the cortex from below are also modeled by a (smaller) circle, here just touching the larger surrounding circle. Any convenient geometry will serve our purposes here: from the simplest (a semicircle) or a general crescent or lune, or as we illustrate in Figure 2(a).

The detailed form of the energy density used is described in Appendix A. In addition, for numerical purposes, one could consider a polygonal (piecewise-linear) approximation to the curve as shown in Figure 2(b), which is explained further in the appendix. For the purposes of infinitesimal symmetries, we use the continuum version.

2.2.   Growth, form and symmetries

We now turn to the approach to the formalism for obtaining the generators and invariants of growth. As mentioned above, ${\bf{u}}\left(s\right) \; \in\mathbb{E}^{2}$ is the configuration of a curve in the plane parametrized by the number $s\in\lbrack0, 1]$. The stretch of the curve is $l\left(s\right) : = \left\vert \frac{d}{ds}{\bf{u}}\left(s\right) \right\vert$. Then the reference metric in question are reference values for $l\left(s\right)$ denoted by a subscript, e.g., $l_{0}\left(s\right)$. Growth is driven by changes with time of this reference metric. The shape is obtained by minimizing energy at each time point. Both the energy and the reference metric can depend on the configurational variables, their derivatives, and even non-locally on their integrals. This approach retains its validity in higher dimensions where we have a metric tensor in the reference state as well as the deformed state. As is well known the strain tensor is just the difference between the pullback to the reference state of metric of the deformed state, and that of the reference metric.

2.2.1.   A toy example

We first describe a toy example that illustrates what we propose, followed by the equations that arise from the model alluded to above and described in Appendix A. Consider variables $x, y$ and a parameter $s$, $0\leq s\leq2\pi$ so that a circle of radius $1$ in the plane is described by

Let us define an energy function $E$

This assumes its minimum along the helical "space curve" whose projection onto the $x-y$ plane is the circle. Now any transformation of the circle to another circle, or indeed closed curve in the plane, can be described by the same range of parameter $s$, so one can compare the positions $x, y$ at the same value of$s$, before and after the transformation. If we introduce the infinitesimal changes

then for example the value of the transformed field $x^{\varepsilon}$ at the original value of the parameter is

and similarly for $y$. Using Eq (2.5) for $x^{\varepsilon}$and the corresponding equation for $y^{\varepsilon}$, the derivative of the energy function at the original circle is

Thus if $Q_{x}$ and $Q_{y}$ are "arbitrary" variations, then the derivative of the energy vanishes only on the circle mentioned. If we compute the second derivative restricted to the circle

we get :

which is the quadratic coefficient of an expansion of the "energy" around the circle. Demanding this be a minimum results in that $x\rightarrow \; x+\varepsilon Q_{x}, y\rightarrow \; y+\varepsilon Q_{y}$ transforms circles to circles provided

Substituting $x = \cos s$, $y = \sin s$, we get

In polar coordinates, $\rho, \theta$, we could write the infinitesimal generator as

and the condition of symmetry becomes

The usual solution offered to describe the symmetry of the circle is the solution

giving the infinitesimal generator of rotations in the plane in the form $y\frac{\partial}{\partial x}-x\frac{\partial}{\partial y}$, or in polar coordinates, $\frac{\partial}{\partial\theta}$. The eigenfunctions of this, namely the trigonometric functions or periodic exponentials, are useful in analyzing functions on the circle. However, any $T$ will describe at least a local symmetry of the circle. For example if we choose

this generates rotations of up to an angle of $\pi$, starting for example at $\theta \; > \; \frac{\pi}{2}$, and up to $\theta \; < \frac{3\pi}{2}$, points on the $y-$axis being "infinities" of the transformation. By a similar token,

gives us a flow or vector field which has singularities in the form of stagnation points instead of infinities, all along the $y-$axis. We believe that such restrictions on the exponentiation of Lie algebras is the norm in biological applications, where the symmetries of shapes are rather local. The same result can be obtained from the symmetric bilinear form in the variations $x\rightarrow x+\varepsilon U_{x}+tQ_{x}$ and similarly for $y$. We consider now "the problem of equivalence", or "growth" as we will understand it, driven by a change of the parameter that describes the particular circle, namely the radius. Denoting

the bilinear form ${\bf{u}}^{\left(1, 1\right) }$ evaluated on the circle is

The first variation is simply to obtain the equation describing a circle, while the second is to obtain the transformation that indicates the growth of the circle of radius $1$ to one of radius $1+t$. Equating to zero the coefficients of the arbitrary variations gives us

which yields

from the second variation. This is the condition for the coordinates $x, y$ of the locus points to remain on a circle when its radius is increased. Again, in polar coordinates, we now have the condition $R = 1$. The "special functions" associated with the operator $\frac{\partial}{\partial\rho}$, namely the real exponentials, are useful in describing radial expansions. We note that $T$ (the coefficient of the vector field along $\theta$) is not determined, and can be arbitrary, being any solution of the homogeneous form of the equation. However, we are quite uninterested in $T$, since any $T\frac{\partial}{\partial\theta}$ is a symmetry common to both the original and the "grown up" circle. Among the equivalence transformations therefore, we must weed out the common symmetries. We can do so if we know these common symmetries, by projecting onto the subspace – at each tangent space – orthogonal to the generators of the common symmetries: the projection of $T\frac{\partial}{\partial\theta}$ is zero along the direction of the vector $\frac{\partial}{\partial\rho}$ at every point $(\rho, \theta)$ of the plane. Further, in general, we will also need to discard the trivial equivalences which, in this two-dimensional case, are of the form $X = 0 = Y$; i.e., where the change in parameter driving growth, the radius of the circle in this example, can be accommodated without transformation of the dependent or independent variables. This does not occur in our particular example, but can in general. Another remark on this example: without further descriptions of the parameter changes, we cannot alter the form of $R$ above. As an example, we cannot describe an upper limit to the radius without putting it in by hand. The change in parameter has been "parachuted" in as a source, being quite unaffected by the disposition of the points of the curve. In general, the transformation of the parameter (the concentration of a chemical for example) will depend on parameters of the shape (the curvature of the cortex for example). Finally, we should remark that although we have, in this example, talked about vector fields in the plane, they need to be defined only on the circle being considered. However, in general they will point "out" of the manifold in question so they can generate growth in the embedding space.

Thus our proposed program in general is: study the equivalences, after barring the common symmetries and trivial equivalences, that map solutions of equations that contain parameters, to solutions of the same equations with altered parameters. The general methodology was apparently introduced by ^[115], Chapter $6$, and developed considerably by Ibragimov and others. For example, Ibragimov, Torrisi and Valenti ^[116] find the equivalence group for a large class of equations. A detailed development of these methods for a general system of differential equations not necessarily arising from a variational principle as in our case, with coefficients that can be arbitrary functions of one or more variables in the problem, may be found in the thesis ^[117]. The thesis also describes the terminology of common symmetries and trivial equivalences. As also indicated in these references, the resulting equivalences will be described by a Lie algebra that will usually generate only a local group.

This conclusion may be summarized in the following diagram:

In the top row, a transformation of parameters changes the equations of equilibrium. The left vertical descent shows a symmetry of the equations of equilibrium prior to a growth step, with $P$ as its generator. The bottom row shows a map from a solution of the prior equilibrium equation to the post-growth equilibrium equation, with $Q$ as the generator. The right vertical descent shows similarly a symmetry generator of the post-growth equation, which is naturally related to the generators of the pre-growth symmetry as well as the that of the transformation of the equation itself.

Below we describe our point of view on growth as a succession of equilibrium states. We will always stay within the infinitesimal, hence a local, look at the shapes. In other words, we will not integrate the algbera to obtain a group. We also feel it is essential in the biological context to deal with Lie groups that may be only local. Special functions that we hope turn out to be useful in describing the shapes generated may be associated to these algebras. When the equations are derived from an extremization principle, this study can be carried out in the context of a second variation. More generally, we could carry out the suggested program in terms of equivalence transformations of the equations themselves. However, in this paper we will confine ourselves to equations arising from a variational principle.

We can now proceed to the general case, which is only notationally but not conceptually more complicated than this toy example.

2.2.2.   Application to the cortical contour

In Appendix A, we describe the detailed form of the energy terms that add up to give the total energy of Eq (2.1). Thus ${\bf{u}}\left(s\right)$ is the location of a point on the 'cortical contour' embedded in the plane. In this section, we assume some familiarity with ^[111] and its notation, and particularly Chapter $5$ in that book. In order to write formulas compactly, we define the derivative of $\mathcal{E}$ at ${\bf{u}}$ evaluated on ${\bf{U}}$

We adopt the usual notation of a linear operator acting on a vector by juxtaposition to denote the action of the derivative on the tangent vector ${\bf{U}}$. Although it is important for consistency to denote ${\bf{u}}$ as a subscript in the derivative operator, we sometimes omit it for typographical ease. Note that

where the angle brackets indicate the integral of a product of a constant function $\not 1$ of the jet space variables with $E^{\prime}{\bf{U}}$. The adjoint of the derivative can be uniquely defined by integration by parts or by

$\left(E^{\prime}\right) ^{\ast}\not 1$ is the Euler Lagrange expression. We will use explicit formulas in terms of differential operators for all these derivatives later. By extension, the second derivative of $E$ is

Clearly $\mathcal{E}^{\prime\prime}$ is symmetric from its definition, and it does not matter which one of ${\bf{U}}, {\bf{Q}}$ go into the first space to the right of $E^{\prime\prime}$. By differentiating (2.24)

Numerical subscripts will always denote the order of a derivative with respect to the independent variables, never components of a vector. Thus:

The use of the letter $J$ to stand generically for the derivatives can be generalized to be a multi-index when several independent variables are involved, as is done systematically in Olver's book. Variations will follow the same convention for derivatives, e.g., if ${\bf{u\rightarrow u}}+\varepsilon{\bf{U}}$, then ${\bf{u}}_{J\;}{\bf{\rightarrow u}}_{J}+\varepsilon{\bf{U}}_{J\;}$. A function of several such differential functions can be written, for example, as

We proceed as we did in the trivial case of the circle above. We define two-parameter transformations ${\bf{u}}\left(s\right) \rightarrow {\bf{u}}\left(s\right) +\varepsilon{\bf{U}}+t{\bf{Q}} \; : = {\bf{u}} ^{2}$ and similarly for any other functions. We use the following notation:

We have deliberately omitted the arguments of the functions and their variations, because this will be the subject of subsequent discussion, and choice of those arguments on which the variations depend will directly affect the obtained results. The notation for the derivatives uses subscripts: ${\bf{u}}_{1}^{\varepsilon}$ means $D{\bf{u}}^{\varepsilon}$ and similarly. For the length density we introduce

A change of the length density is what drives growth in this case. Finally, we introduce the further special notation in this section:

That is, all the relevant arguments in these differential functions take on their values at non-zero $t$. Recall that

and that the total derivative operator is here:

As for example ^[111], Chapter 5 indicates, if one allows the generators to depend on derivatives, then it is sufficient to transform the dependent variables and leave the independent ones alone. For the sake of completeness, we recount the argument briefly. Consider the solution near a particular $s$

and apply the operator $S\partial_{s}+{\bf{V}}\cdot\partial_{{\bf{u}}}$ to obtain

as a condition of symmetry. This is the same as applying ${\bf{Q}} \cdot\partial_{{\bf{u}}}$ on solutions, where

with ${\bf{u}}_{1}\leftarrow{\bf{f}}^{\prime}\left(s\right)$. Since all our considerations are local, the representation (2.35) is general enough. Another argument is that as long as the domains of the independent variable remain identical before and after the transformation, there is no loss of generality in comparing the transformed dependent variables with those of the untransformed at the same value of the independent variable. Then

where

The superscripts $t$ remind us that the corresponding expressions are being evaluated at the extended functions ${\bf{u}}^{t}$, $l^{t}$; i.e., $t$ has not yet been set to zero. Continuing with the expansion,

where the neglected terms in the last equation are neither $\mathcal{O}\left(t\right)$ nor $\mathcal{O}\left(t\varepsilon\right)$. We set the $\mathcal{O}\left(\varepsilon\right)$ variation to zero, so that necessarily ${\bf{\Delta}}E\left({\bf{u}}, l\right) = 0$.

Remark 1. The first variation in ${\bf{Q}}$ (i.e. $\mathcal{O}\left(t\right)$ term) does not vanish on $ {\bf{\Delta}}E = 0. $

Continuing, we compute the terms in the second variation above separately. First,

Explicitly

We will now assume that $l$ is given as a function of $s$ alone, and has no explicit dependence on${\bf{u}}_{J}$. Then, since $\left[ \hat {v}_{{\bf{Q}}}, D\right] = 0$ we have, denoting the Hessian of $E$ in the ${\bf{u}}$ variables by $H_{IJ}^{{\bf{u}}}$, that

Similarly

When "dotted" into the variation ${\bf{U}}$, $\hat{v}_{{\bf{Q}}}$ and $\hat{v}_{L}$ comprise the important part of the second variation. All of the above commutations will fail if the reference length functions $l$ depend on differential or other functions of the dependent variables ${\bf{u}}$ (or vice versa). In that case the first line of the expressions for $\hat{v}_{{\bf{Q}}}\left({\bf{\Delta}}E\right)$ must be retained. We proceed with the simplifying assumptions, where $l$ is given as a function of $s$ alone.

Now let us return to the boundary terms. The second variation is:

where the sum truncates according to the order of the differential function it operates on. Denoting $\hat{v}: = \hat{v}_{{\bf{Q}}}+\hat{v}_{L}$, and noting that ${\bf{U}}$ vanishes at either boundary,

Again, by assuming that ${\bf{U}}$ is not a differential function in the ${\bf{u}}_{J} \; ^{\prime}$s, $\hat{v}$ annihilates ${\bf{U}}, D{\bf{U}}$ so

The coefficient of ${\bf{U}}_{1}$ in $\mathcal{B}_{1}$, namely $\frac{\partial E}{\partial{\bf{u}}_{2}}$ must vanish at the end points, from the fact that the first derivative of the energy evaluated on ${\bf{U}}$ must vanish at the extremal. Further simplification of $\mathcal{B}_{2}$ will depend on the special forms for the energy density. When we examine the second variations from the point of view of the analysis of equivalence, we do not find that the boundary terms will vanish, owing to the more general forms of the variations that we need. These vector fields which minimize the second variation are called Jacobi fields in the calculus of variations and Riemannian geomety. Jacobi fields are usually discussed in the context of a metric which is a quadratic form in the "velocities" (the usual Riemannian space), and in the context of how "fast" geodesics separate along a path. In the present case, we are talking about the separation of an extremum for one energy function from an extremum starting at the same point but for another energy (the one with a change of reference metric). For this reason, the first variation of the energy on this field does not vanish. Then ${\bf{Q}}$ too must satisfy the Euler-Lagrange equations if ${\bf{u}}$ does and if the transformed curve is also to be a solution for the new rest length ($l+tM$). The infinitesimal version then means

Since ${\bf{\Delta\cdot U = }}\mathcal{E}_{{\bf{u}}}^{\prime}\left[ {\bf{U}}\right]$, this means

This gives us integro-differential equations involving the fields $S, {\bf{u}}, L$ and their derivatives with respect to $s$, which we may attempt to solve to give us the "vector fields that generate growth". However, the Eqs (2.54) or (2.55) are enormously complicated. The gradients and the Hessians of the energy densities appear in the equations, as is seen from Eq (2.1). These are evaluated in Appendix A, but the total derivatives $D$ still need to be expanded. As an example,

Such expansions are feasible using computer algebra packages, but the resulting equations that need to be solved to obtain the generators are very complicated. The framework nevertheless may prove useful in the future, but for now we turn to a rather more immediate possibility for exploration.

3.   Lie algebras from data

Fortunately, there are sources other than the mathematical models expressed either as differential equations or as minimization problems to obtain approximations to the vector fields that will serve well for numerical purposes. We shall illustrate how we may use them in the following section 4. In this section we give several examples, both from our own work and those of others, to illustrate the wealth of data available to construct growth generators. First we should mention that all the models that people have developed (see section 1) may in fact also be considered as potential sources of the empirical vector fields in question, since one may use them to examine the movement of corresponding points of the brain during growth. We discuss a specific example below. We shall order the examples in decreasing dimensionality, which corresponds roughly to increasing simplicity. In this section, we illustrate some of the data available. Many details about the datasets including where they may be accessed are contained in Appendix B and the links there provided to Supplementary Information.

3.1.   Three-dimensional vector fields in a volume

We obtained the original $16$ datasets from a research project described in ^[15]. We then processed the data to obtain vectors that point from a position at one age to the corresponding position at the succeeding age as outlined in Appendix B. The process involved both global (affine) and local (so-called elastic morphing) transformations to co-register into a common frame any pair of volume data to be able to obtain correspondence between the voxels of the pair. Vectors that characterize the growth patterns are available in the NIHto682dispVector folder. Further interpolation results in the data in the folder $NIH_GrowthVecDataInterpolated$ and these datasets are available for exploration. Figure 3 shows a selection from the vast amount of data. We selected a particular plane and displayed vectors that lead from one time point to the immediately succeeding one at the voxels in that plane, at four different time points separated roughly by one year. The displays are color-coded line-fields as opposed to vector fields, i.e., the magnitude of each component (with sign ignored) of the growth vector is displayed as an intensity of a corresponding R, G, B color.

We also created volume renderings of the magnitudes of the vectors. In Figure 4, we show select volume renderings at roughly one year apart. The figure shows still frames from an animation which may be viewed from the repository described in Appendix B. The actual dates selected were those of the datasets designated $0682, 1090, 1552,$and $2058$. Precise dates for these are given in the Appendix.

The data above was that of a post-term child, and it is known that there is essentially no gyrification (increased folding of the cerebral cortex) beyond the age of about one year (see references mentioned in the Introduction). It is therefore of considerable interest to examine the brain growth vectors in pre-term infants.

3.2.   Three-dimensional vector fields on the cortical surface

A superb set of MR images, enabled with point to point correspondence has been made available publicly ^[13]. These were taken of generally pre-term infants at birth at PMA (post-menstrual ages) of 27, 31, 34, and 37 weeks, respectively. We have constructed the coordinates of the corresponding points based upon the information provided by the researchers, and these are available as Excel spreadsheets at the location mentioned in Appendix B, which also contains more useful information on these data. A plot of the vector fields is shown in Figure 5. These are plots of the vectors connecting the corresponding points taken from the spreadsheets provided. We were provided with two different methods of registration between data at different points, and for the purposes of the display here, we chose the forward registration as mentioned in Appendix B. With this choice, Figure 5(a) is a plot of vectors that lead from the cortical surface at week $27$ to the surface at week $31$, and analogously for the other figures. The data provided contained 10,242 points on each surface, and we have plotted the vectors at every $15^{\rm{th}}$ point. The pattern of growth which reverses itself to return to point to the rostral region of the brain was noted in the earlier work cited ^[15] as a "rostro-caudal wave". We note that sub-Figures (a)–(c) show the vectors between consecutive growth stages imaged, while Figure 5(d) shows the vectors from week $31$ to the very last (at-term) image at week $37$.

From the spreadsheets shown in the link referenced in Appendix B, other vectors, e.g., arising from a "backwards" registration may also be obtained and examined (see the Appendix for details). The above plots were created from the Excel files mentioned by use of the Vector command in Mathematica 12.1.

3.3.   Two-dimensional vector fields on a curve

The Ph.D. thesis of Sarah Kim ^[91] develops a model for the cortical curves that arise from sectioning the brain. Actually, these are 'thickened' curves defined by quadrilateral elements. The model she provides then moves the vertices of these quadrilaterals, and one may see from Figure 4.3 of her thesis that the details within her numerical simulations allow one to construct the vector fields of growth. Indeed they can be defined throughout the thickness of the cortical layer by interpolation between the outer and inner vertices of the quadrilaterals. Both her model and our previous work ^[114] construct only these boundary curves, though thickened in her case. However, her model is more general than ours because she allows for two dimensional interactions. In our previous work, we invented a penalty function or energy to prevent self-intersection of the cortical curve and allow for white matter fibers to pass between points on the cortical curve. This allowed us to stay within the mechanics of curves in the plane, without having to take into account any interaction between points on the curve and those outside. We took advantage of special geometry of planar curves, which would not be possible in a higher dimension.

There is another potential source of such fields, namely extensive data collected on brain sections. The monumental six-volume set authored by Shirley Bayer and Joseph Altman ^[118] contains pictures of extensive collections of such sections taken from the first trimester to term. In Figure 6, we show a few of these sections which also serve to highlight the problems in using such un-registered data collected historically. The sections shown are from what is called the Yakovlev collection ^[119] at gestational weeks $29$, $31$, $34$, and $37$, not dissimilar to the data mentioned above. Moreover, we have selected the sections that most closely correspond at the different weeks.

Despite this, we can easily see that it is difficult to obtain a detailed point-to-point correspondence and be confident of not contaminating the vector fields with artefacts due to poor correspondence and registration. Our own early work also readily provides empirical vector fields by following marked points in growing curves, but that model was more purposed to look at the basic determinants of growth rather than detailed correspondence with cortical curves.

As mentioned, this section was primarily to indicate the wealth of sources available, and continuing to be available for vector fields of both two dimensions in the plane or on curves, and of three dimensions in a volume or on a surface.

4.   From generators to special functions

We now return to our main theme and speculation, namely that Lie algebras and special functions associated with them may be useful for describing the mathematical development of biological pattern formation and growth. Let us henceforth specialize to the simplest case, namely the growth of a cortical curve, or the cortical boundary of a plane section as described above. We may locally describe a growth transformation as taking a material point at $(x, y)$ and moving it to $(x+ \; \varepsilon, \; y+\varepsilon q\left(x\right))$. We then call $q\left(x\right)$ the generator of the growth pattern, and we have supplied many examples in diverse dimensions. Now, we can construct the Lie algebra generated by $\partial/\partial x$ which we shall henceforth call $D$, and multiplication by $q\left(x\right)$. Sometimes, it is convenient to introduce a parameter $t$ and call the generators $tD$ and $q\left(x\right)$. If $q\left(x\right)$ is a polynomial this generates a finite dimensional nilpotent Lie algebra which includes multiplication by $q\left(x\right)$ and its successive derivatives. The canonical prototype is the Heisenberg algebra in which $q\left(x\right) = x$. Its $N-$step nilpotent generalization is $q\left(x\right) = x^{n}$ or $x^{n}/n!$. If $q\left(x\right)$ is smooth and infinitely differentiable, then an infinite dimensional Lie algebra, though with only two generators, will result. The following developments can easily accommodate either case.

There are many roads available to take us to special functions. To construct some of these, and display them as graphs, we shall choose $q\left(x\right)$ to be the bump function so labeled in Figure 7(a). We also take the range of $x$ to be a finite interval, $[-1, 1]$ without loss of generality. This is not strictly necessary, but allowing for infinite ranges will mean that we have to pay attention to the integrability, and we choose to ignore that for the moment. Moreover, in the examples we have illustrated, the skull is a compact bounded region, and so indeed the ranges will be finite. Now, we notice that $D$ and $q\left(x\right) -D$ are operators adjoint with respect to each other in the space of functions which vanish at the end points, and where the inner product has the weight $w\left(x\right) \sim\exp\left(-Q\left(x\right) \right)$, $Q^{\prime}\left(x\right) = q\left(x\right)$. The function $Q\left(x\right)$ and the weight $w\left(x\right)$ after multiplication by a factor ensuring normalization, are also displayed in Figure 7(a). Immediately then, the "growth generator" gives rise to an orthogonal polynomial sequence with the noted weight function. We construct these by standard procedures, and no particular difficulties arise doing so numerically when $q$ and hence $Q$ and $w$ are only given numerically or 'empirically' as we have called it. The first five orthogonal polynomials for the given $q$ are shown in Figure 7(b).

We may develop other special functions. We note that the product $D \; \left(q(x)-D\right)$ is a self adjoint operator and hence its eigenfunctions in the space spanned by functions vanishing at the end points may also be considered "special". The first five of these are plotted in (c) of the figure. The pair $D, q(x)-D$ of operators generate the same Lie algebra as do the original pair $D, q(x)$ or other pairs such as $x, q(D)+x$. Finally we may consider the Hamiltonian operator with a potential well $-q\left(x\right)$. The eigenfunctions of $-D^{2}-q\left(x\right)$ are shown in Figure 7(d). All these functions plotted may be considered "special functions" associated with the generators of the Lie algebra in question. Another set of operators whose eigenfunctions are sometimes designated as special functions associated to the Lie algebra are the quadratic Casimir operators. However, in our case, these are not useful because they do not involve the differentiation operator. There do exist such operators in an isomorphic Lie algebra which are quadratic in the Lie algebra elements but these are of high order in $D$. See Appendix C for further discussion about the Casimir operators.

We have not yet made any use of the fact that these generators generate an algebra, and that other Lie algebra operators are thereby obtained. As stated, even in the infinite dimensional case, we may consider the whole class of elements of the Lie algebra as generalizations of the standard $N-$step nilpotent algebras. The relation between Lie groups and algebras and special functions is a venerable subject, but standard texts usually focus on semisimple Lie groups. The papers of Feinsilver and Schott offer methods seemingly applicable to any Lie algebra, and moreover focus on the algebras, rather than the groups. In particular their methods apply to nilpotent and solvable Lie algebras, and algebras generated as above are clearly either nilpotent or obvious generalizations of such algebras. We have consulted the papers ^[120,121], and ^[122] and we should also mention their three-volume tome ^{[123,124,125]} summarizing several decades of their work in the field. For the reader's convenience we have provided a short exposition of Feinsilver and Schott's methods in Appendix C. The Lie algebra representations on suitable spaces of functions of $x$ generated in the special case mentioned consists of the multiplications by $q\left(x\right), q^{\prime}\left(x\right), .....q^{(n)}\left(x\right), ....$ and the differentiation operator $D$. When $q$ is infinitely differentiable, the Lie algebra may be infinite-dimensional. The basic machinery that Feinsilver and Schott use to generate "special" functions and polynomials is to write an element of the Lie group in two different ways. For any suitable finite dimensional Lie algebra with elements $X_{i}$, $i = 1, 2, ..., n$, we exponentiate to write an element of the Lie group so obtained in two different ways

This equality gives us the parameters $a$ in terms of the $\alpha^{\prime}$s, called coordinates of the first and second kind. The details of how these may be obtained in terms of the adjoint representation of the Lie algebra is given in Appendix C. Let us take the case of the $4-$step nilpotent Lie algebra. Let the non-zero commutators in the Lie algebra be written as

The generators are $X_{5}$ and $X_{1}$; for example $X_{5} = tD$ and $X_{1} = x^{3}/3!$ generates such an algebra and $X_{4}$ is then the center of the algebra. We find in this case

The general pattern is clear from this set of equalities, and a proof that does not depend on working out the transformation of coordinates this way is indicated below. Further, it is clear from the derivation that the relationship does not depend on the specific representation chosen and that any morphism of the Lie algebra will yield the same relationships. The general procedure for obtaining special functions suggested by the work of Feinsilver and Schott, specialized to the specific representation of the Lie algebra above, is this. Apply both sides of the Eq (4.1) to some suitable function, for example the constant function $1$. There results on the right hand side the expression

since $\exp\left(\alpha_{5}X_{5}\right) 1 = \exp\left(\alpha_{5}D\right) 1 = 1.$ It is assumed that we have replaced $a_{i}^{\prime}$s by their corresponding expression in terms of the $\alpha^{\prime}$s as well as the $X_{i}^{\prime}$s in terms of the multiplicative functions. Then we may examine the moments for any multi-index $J$ given by

As implied,

Let us take a particular example. This example when restricted to the case $n = 3$, which is the Heisenberg algebra, will yield the Hermite polynomials as the moments. We will take this example for any $n$ including $n = \infty$. Let us focus on the generators only, setting $\alpha_{1} = \alpha$ and the $\alpha$ for the differentiation operator to be $\beta$. The other coordinates are set to zero for the moment. All of the $X_{i}$ in Eq (4.4), now with the upper limit at infinity, commute, and we obtain

where we remind ourselves that $Q^{\prime}\left(x\right) = q\left(x\right)$ and finally setting $\beta = -\alpha$ we obtain the generating function

Setting $Q\left(x\right) = x^{2}/2$ generates the Hermite polynomials, as mentioned. Since

we have

Thus $\Phi$ generates the moments

This generating function and its attendant moments may also be considered as functions of interest allied to the Lie algebra. In fact, we can generalize this approach without setting the other parameters to zero. Consider an element of the algebra

where the equations define the functions $\psi, \Psi$. Then

and we have generating functions for other potentially interesting functions associated with the algebra. We now use this to obtain the general relationship between the two kinds of coordinates, which we surmised by inspection of a lower order case in Eq (4.3). By Taylor expansion in $t$, we find

Setting $t = 1$ reproduces the special cases we considered before.

We have thus shown a variety of "special" functions associated with the generators of the Lie algebra. Although the plots were specific to a model bump function, it is clear that any empirical data of the proper sort can be used to obtain a generator which in general will be a vector field in space, and will generate a Lie algebra in the usual way. The supplementary file referred to above contains several other morphisms of the Lie algebra which may lead to interesting functions and polynomials.

5.   Conclusions

This paper has consisted of three parts. In the first part, we model early brain growth as smooth deformations (diffeomorphisms) of a spatial region containing the brain which gives rise to cortical folds. The rate of change of these diffeomorphisms are what we have called the vector fields of growth. We applied the approach of symmetries of equations, in particular those derived from a variational principle, to some simple elastic models of growth to suggest equations that may describe the generators and invariants of growth, specifically of the brain and cortex. The equations for even a simplified one dimensional model were very complicated. We then showed how vector fields that are the growth generators may be obtained from a variety and increasing number of sources, both from computer models and from imaging. We hypothesize that these vector fields are finitely generated. In the third part we showed new special functions that can be obtained from such generators. We emphasize that the methods described may be extended to higher dimensions and more complicated generators than the one used for illustration in the previous section. The hypothesis that a finitely generated algebra is useful and sufficient may be tested by examining the vector fields of growth at different times and confirming this space of vector fields is approximately finitely generated. This development from empirical data to Lie algebras with associated special functions is the principal idea in the paper. It is our hope that developments of the work presented here will be useful in creating new approaches to computational modeling and prediction of brain growth and to improve qualitative understanding of the physiological drivers of growth and form.

Acknowledgments

The author thanks a number of people for help. The data on growth vectors listed under the description Three dimensional vector fields in a volume in Appendix B was created with the help of "elastic" registration programs kindly made available by Rainer Lachner of BrainLab AG. Kara Garcia of the Indiana University School of Medicine kindly supplied a new set of surface files with optimized and fewer points than those available from the published paper as Appendix B describes. Also as mentioned there, Shirley Bayer of Ocala, Florida made available her slides from the Yakovlev collection which we have referenced. Professor Philip Feinsilver of the School of Mathematics in Southern Illinois University explained his techniques which are summarized in the Appendix C. Two anonymous reviewers pointed out some errors in the first submission. Their summaries of this work have been used in the Conclusions section as well. The author is responsible for all the work reported in the paper and any errors therein, unless acknowledged otherwise.

Conflict of interest

The author declares that he has no conflict of interest with regards to this work.

Appendix

A.   One-dimensional cortical contour model

We briefly describe the energy density terms in the cortical contour model described in ^[114]. ${\bf{u(}}s{\bf{)}}$ is a 2D vector with components $(u_{1}(s), u_{2}(s))$ which denote the $x-$and $y-$ components of the position of the curve in the plane in parametric form, parametrized by $s\in\lbrack0, 1]$. The prime denotes differentiation with respect to $s$. ${\bf{K}}$ denotes the curvature vector (i.e., the acceleration of ${\bf{u}}$ with respect to a variable representing arc-length parametrization). $l$ in the White matter entry refers to length along the curve, while $\vert \; \vert$ denotes the usual Euclidean distance in the plane. $\Theta$ is the usual Heaviside function, i.e, 1 when its argument is positive and $0$ otherwise. The point denoted $c$ is an approximate "center" of the curve, or that of its approximately semicircular ceiling. This point could be taken as the origin of coordinates.

In the discrete model there are a set of nodes numbered $0$ to $N+1$. The various energy terms must be summed over the number of variables. For example, the stretching term must be summed from $i = 1$ to $i = N$. Referring to Figure 2(b), the lengths $l$ are those of the polygon sides, and the angles are those between adjacent sides. There are many equalities that connect the various variables, all from elementary plane geometry.

The total energy is an integral over $s$ in the interval $[0, 1]$ of the sum of the terms listed above, and the actual curve is obtained as follows. First an initial curve is prescribed with total length equal to $\int_{0}^{1} dsl_{0}\left(s\right)$. Then, in an actual placement of the curve in the plane, the energy density (energy per unit parameter $s$) is the sum of the terms listed above. The first term is the energy of stretching with a non-zero rest length (per unit $s)$, and is just Hooke's law. The second is an energy of bending, again with a stress-free state with curvature ${\bf{K}}_{0}$. The multipliers which are the signs of the curvatures are needed to make the energy density invariant with respect to rotations. As can be seen, there are non-local terms, or double integrals in the expression of the energy, despite this being a model using mechanics by contact. This is because we can approximate the effects of the white matter "cables" without modeling these explicitly. As is often the case, when we eliminate a set of variables, especially ones that mediate interactions, a local theory of action by contact becomes one of "action at a distance". We didn't actually derive this nonlocal potential from a local theory, but simply postulated it. See ^[114] for details.

A.1.   Gradients of the energy density

Below we list the (non-zero) gradients of the energy density with respect to the dependent variables and their differentials, as needed to write down the Euler-Lagrange conditions for the extremals. The self-avoidance term mentioned in Table 1, which would be $E_5$, is omitted since it is not needed in the model when the terms $E_6$ and $E_7$ are included.

$E_{1}: = \; [|D{\bf{u}}|-l_{0}]^{2}/l_{0}^{2}$

$\frac{\partial E_{1}}{\partial\left(D{\bf{u}}\right) } = \frac{2}{l_{0}} \; \left(1-\frac{l_{0}}{|D{\bf{u}}|}\right) D{\bf{u}}$

$E_{2}: = \int_{0}^{1}ds \; (\sigma|{\bf{K}}|-\sigma_{0}|{\bf{K}}_{0} |)^{2} = \int_{0}^{1}ds\left(K\mp K_{0}\right) ^{2}$ according to whether the sign (well defined for planar curves) of the two curvatures are the same or not.

$\frac{\partial E_{2}}{\partial{\bf{u}}_{1}} = 2\sigma\left(K\mp K_{0}\right) \left(-\frac{\mathcal{J}{\bf{u}}_{2}}{\left\vert {\bf{u}}_{1}\right\vert ^{3}}-3\frac{|{\bf{K}}|}{\left\vert {\bf{u}} _{1}\right\vert ^{2}}{\bf{u}}_{1}\right)$ where $\mathcal{J}\left(u_{x}, u_{y}\right) : = \left(-u_{y}, u_{x}\right)$

$\frac{\partial E_{2}}{\partial{\bf{u}}_{2}} = 2\sigma\left(K\mp K_{0}\right) \frac{\mathcal{J}{\bf{u}}_{1}}{\left\vert {\bf{u}} _{1}\right\vert ^{3}}$

$E_{3}: = \int_{0}^{1}ds (d-u_{y})^{2}\times\Theta(d-u_{y})$

$\frac{\partial E_{3}}{\partial u_{y}} = 2(u_{y}-d)$ is the only non-zero component of the gradient of $E_{3}$.

$E_{4}: = \int_{0}^{1}ds \; (|{\bf{u}}-{\bf{c}}|-r_{0})^{2}\times \Theta(|{\bf{u}}-{\bf{c}}|-r_{0})$

$\frac{\partial E_{4}}{\partial{\bf{u}}} = 2\frac{\left\vert {\bf{u}} -{\bf{c}}\right\vert -r_{0}}{\left\vert {\bf{u}}-{\bf{c}}\right\vert }\left({\bf{u}}-{\bf{c}}\right)$

$E_{6}: = \int_{0}^{1}ds_{1}\frac{\left(\int_{s_{1}}^{s}ds_{2}\left\vert D_{s_{2}}{\bf{u}}\right\vert \right) ^{p}}{\left(\left\vert {\bf{u}}\left(s\right) -{\bf{u}}\left(s_{1}\right) \right\vert +\epsilon\right) ^{p}}$

$\frac{\partial E_{6}}{\partial{\bf{u}}} = - \; {\bf{1}}p\int_{0}^{1} ds_{1}\left[ \frac{\left(\int_{s_{1}}^{s}ds_{2}\left\vert D_{s_{2} }{\bf{u}}\right\vert \right) ^{p}}{\left(\left\vert {\bf{u}}\left(s\right) -{\bf{u}}\left(s_{1}\right) \right\vert +\epsilon\right) ^{p+1}\left\vert {\bf{u}}\left(s\right) -{\bf{u}}\left(s_{1}\right) \right\vert }\right] +$

$\frac{\partial E_{6}}{\partial{\bf{u}}_{1}} = p(p-1)\frac{{\bf{u}} _{1}{\bf{u}}_{1}}{\left\vert {\bf{u}}_{1}\right\vert ^{2}}\int_{0} ^{1}ds_{1}\left[ \frac{\left(\int_{s_{1}}^{s}ds_{2}\left\vert D_{s_{2} }{\bf{u}}\right\vert \right) ^{p-2}}{\left(\left\vert {\bf{u}}\left(s\right) -{\bf{u}}\left(s_{1}\right) \right\vert +\epsilon\right) ^{p} }\right] +$

$E_{7}: = \int_{0}^{1}ds [\frac{d}{ds}|{\bf{u}}-{\bf{c}}|]^{2n}$

$\begin{array} [c]{c} \frac{\partial E_{7}}{\partial{\bf{u}}} = n\frac{\left(\left[ {\bf{u}}-{\bf{c}}\right] \cdot{\bf{u}}_{1}\right) ^{2n-1}} {|{\bf{u}}-{\bf{c}}|^{2n}}\left[ D{\bf{u-}}\frac{\left(\left[ {\bf{u}}-{\bf{c}}\right] \cdot D{\bf{u}}\right) }{|{\bf{u}} -{\bf{c}}|^{2}}\left({\bf{u}}-{\bf{c}}\right) \right] \end{array}$

A.2.   Hessian of the energy density

We use what is sometimes called "dyadic" notation below. For example

is a matrix whose components are $\left({\bf{u}}_{2}\right) _{i}\left({\bf{u}}_{2}\right) _{j}, i, j = x, y$. Further, $\mathcal{I}$ is the identity matrix.

$E_{1}: = [|D{\bf{u}}|-l_{0}]^{2}/l_{0}^{2} \frac{\partial^{2}E_{1}}{\partial\left(D{\bf{u}}\right) \partial\left(D{\bf{u}}\right) } = \frac{2}{l_{0}} \left[ \left(1-\frac{l_{0} }{|D{\bf{u}}|}\right) \mathcal{I}+\frac{l_{0}}{|D{\bf{u}}|^{3}}\left(D{\bf{u}}\right) \left(D{\bf{u}}\right) \right]$

$E_{2}: = \int_{0}^{1}ds \; (\sigma|{\bf{K}}|-\sigma_{0}|{\bf{K}}_{0}|)^{2}$

$\frac{\partial^{2}E_{2}}{\partial{\bf{u}}_{1}\partial{\bf{u}}_{1} } = 2\sigma\left(-\frac{\mathcal{J}{\bf{u}}_{2}}{\left\vert {\bf{u}} _{1}\right\vert ^{3}}-3\frac{|{\bf{K}}|}{\left\vert {\bf{u}} _{1}\right\vert ^{2}}{\bf{u}}_{1}\right) \left(-\frac{\mathcal{J} {\bf{u}}_{2}}{\left\vert {\bf{u}}_{1}\right\vert ^{3}}-3\frac {|{\bf{K}}|}{\left\vert {\bf{u}}_{1}\right\vert ^{2}}{\bf{u}} _{1}\right)$

$+2\sigma\left(K\mp K_{0}\right) \left(-\frac{3}{\left\vert {\bf{u}}_{1}\right\vert ^{5}}{\bf{u}}_{1}\mathcal{J}{\bf{u}}_{2} +\frac{6K}{\left\vert {\bf{u}}_{1}\right\vert ^{4}}{\bf{u}}_{1} {\bf{u}}_{1}-\frac{3K}{\left\vert {\bf{u}}_{1}\right\vert ^{2}} \mathcal{I}\right)$

$\frac{\partial^{2}E_{2}}{\partial{\bf{u}}_{1}\partial{\bf{u}}_{2} } = 2\sigma \ \frac{\mathcal{J}{\bf{u}}_{1}}{\left\vert {\bf{u}} _{1}\right\vert ^{3}}\left(-\frac{\mathcal{J}{\bf{u}}_{2}}{\left\vert {\bf{u}}_{1}\right\vert ^{3}}-3\frac{|{\bf{K}}|}{\left\vert {\bf{u}} _{1}\right\vert ^{2}}{\bf{u}}_{1}\right) -\sigma\frac{6}{\left\vert {\bf{u}}_{1}\right\vert ^{5}}\left(K\mp K_{0}\right) {\bf{u}} _{1}\mathcal{J}{\bf{u}}_{1}$

$\frac{\partial^{2}E_{2}}{\partial{\bf{u}}_{2}\partial{\bf{u}}_{2} } = {\bf{0}}$

$E_{3}: = \int_{0}^{1}ds (d-u_{y})^{2}\times\Theta(d-u_{y})$

$\frac{\partial^{2}E_{3}}{\partial u_{y}^{2}} = 2$ is the only non-zero component of the Hessian of $E_{3}$.

$E_{4}: = \int_{0}^{1}ds (|u-c|-r_{0})^{2}\times\Theta(|u-c|-r_{0})$

$\frac{\partial^{2}E_{4}}{\partial{\bf{u}}\partial{\bf{u}}} = 2\left[ \frac{\left(1-\left\vert {\bf{u}}-{\bf{c}}\right\vert \left[ \left\vert {\bf{u}}-{\bf{c}}\right\vert -r_{0}\right] \right) }{\left\vert {\bf{u}}-{\bf{c}}\right\vert ^{2}}\left({\bf{u}}-{\bf{c}}\right) \left({\bf{u}}-{\bf{c}}\right) +\frac{\left\vert {\bf{u}} -{\bf{c}}\right\vert -r_{0}}{\left\vert {\bf{u}}-{\bf{c}}\right\vert }{\bf{1}}\right]$

$E_{6}: = \int_{0}^{1}ds_{1}\frac{\left(\int_{s_{1}}^{s}ds_{2}\left\vert D_{s_{2}}{\bf{u}}\right\vert \right) ^{p}}{\left(\left\vert {\bf{u}}\left(s\right) -{\bf{u}}\left(s_{1}\right) \right\vert +\epsilon\right) ^{p}}$

$\frac{\partial^{2}E_{6}}{\partial u_{i}^{\prime}\partial u_{j}} = -p^{2} \frac{D{\bf{u}}}{\left\vert D{\bf{u}}\right\vert }\int_{0}^{1} ds_{1}\left[ \frac{\left(\int_{s_{1}}^{s}ds_{2}\left\vert D_{s_{2} }{\bf{u}}\right\vert \right) ^{p}}{\left(\left\vert {\bf{u}}\left(s\right) -{\bf{u}}\left(s_{1}\right) \right\vert +\epsilon\right) ^{p+1}}\frac{{\bf{u}}-{\bf{u}}\left(s_{1}\right) }{\left\vert {\bf{u}}-{\bf{u}}\left(s_{1}\right) \right\vert }\right]$

$E_{7}: = \int_{0}^{1}ds [\frac{d}{ds}|{\bf{u}}-{\bf{c}}|]^{2n}$

$\frac{\partial^{2}E_{7}}{\partial{\bf{u}}\partial{\bf{u}}} = \begin{array} [c]{c} n\left(2n-1\right) \frac{\left(\left[ {\bf{u}}-{\bf{c}}\right] \cdot D{\bf{u}}\right) ^{2n-2}}{|{\bf{u}}-{\bf{c}}|^{2n}}\left(D{\bf{u}}\right) \left(D{\bf{u}}\right) \\ -2n^{2}\frac{\left(\left[ {\bf{u}}-{\bf{c}}\right] \cdot D{\bf{u}}\right) ^{2n-1}}{|{\bf{u}}-{\bf{c}}|^{2n+2}}\left\{ \left({\bf{u}}-{\bf{c}}\right) \left(D{\bf{u}}\right) {\bf{+}}\left(D{\bf{u}}\right) \left({\bf{u}}-{\bf{c}}\right) \right\} \\ +2n\left(n+1\right) \frac{\left(\left[ {\bf{u}}-{\bf{c}}\right] \cdot D{\bf{u}}\right) ^{2n}}{|{\bf{u}}-{\bf{c}}|^{2n+4}}\left({\bf{u}}-{\bf{c}}\right) \left({\bf{u}}-{\bf{c}}\right) -n\frac{\left(\left[ {\bf{u}}-{\bf{c}}\right] \cdot D{\bf{u}} \right) ^{2n}}{|{\bf{u}}-{\bf{c}}|^{2n+2}}{\bf{1}} \end{array}$

$\frac{\partial^{2}E_{7}}{\partial\left(D{\bf{u}}\right) \partial\left(D{\bf{u}}\right) } = n\left(2n-1\right) \frac{\left(\left[ {\bf{u}}-{\bf{c}}\right] \cdot D{\bf{u}}\right) ^{2n-2}} {|{\bf{u}}-{\bf{c}}|^{2n}}\left({\bf{u}}-{\bf{c}}\right) \left({\bf{u}}-{\bf{c}}\right)$.

These Hessians enter into the expressions in section 2.2.2. In turn as mentioned there, the total derivatives need to be expanded and finally equations result for the generators $Q$, which if solved, give us the analog of the Lie-symmetry approach to growth and form from mathematical equations thereof.

B.   Supplementary information

We describe here the datasets used and created for our illustrations in section 3. We use the same terminology to aid cross reference. All of the folders referred to by name here are available at https://tinyurl.com/RaghavanSuppInfo at the time of this writing.

Three dimensional vector fields in a volume (NIH data ^[15] used for Figures 3 and 4).

We used 16 MRI datasets obtained with permission from the NIH in-vivo NMR Center, of a child aged 3 years at the starting date of 23 January 1996, and ending on 7 November 2000 (see reference for details and how to obtain these datasets). The data has a starting date of Aug 1996 and the last of the pair is Nov 2000. The specific dates at which these images were obtained are:

1. 0682 Study 2189, Series 2 23 Jan 1996 11:59:28 (3 years old)

2. 810 Study 3458, Series 3 20 Aug 1996 17:47:27

3. 871 Study 3981, Series 2 12 Nov 1996 13:06:09

4. 1019 Study 4911, Series 2 22 Apr 1997 13:13:29

5. 1090 Study 5348, Series 3 01 Jul 1997 11:39:47

6. 1169 Study 5796, Series 3 09 Sep 1997 09:29:19

7. 1251 Study 6513, Series 3 30 Dec 1997 06:50:53

8. 1387 Study 7680, Series 3 09 Jun 1998 13:21:27

9. 1552 Study 8877, Series 3 22 Dec 1998 14:02:36

10. 1662 Study 9468, Series 4 30 Mar 1999 12:47:58

11. 1740 Study 9962, Series 3 16 Jun 1999 13:41:28

12. 1778 Study 10170, Series 2 20 Jul 1999 17:06:22

13. 1930 Study 10983, Series 3 23 Nov 1999 10:15:00

14. 2058 Study 11869, Series 4 28 Mar 2000 13:56:20

15. 2130 Study 12314, Series 3 15 Jun 2000 16:52:34

16. 2284 Study 13116, Series 3 07 Nov 2000 13:11:06 (8 years old)

The datasets in our repository pertaining to these data are in the folders (A) NIHto682dispVector, (B) NIH_GrowthVecDataInterpolated, and (C) NIH_GrowthVecAnimationGIFS.

The original datasets from the NIH were all $256\times256\times124$ volume datasets with pixel size $0.9375\operatorname*{mm}$ and slice separation $1.5\operatorname*{mm}$. These latter numbers are irrelevant for the viewing of the vector data and the displays. The folder A contains 15 volumes, all co-registered to the first "0682" dataset, and consist of the $(x, y, z)$ components of the growth vectors in a $256\times256\times124$ volume. Each component is a floating point real ($32-$bit, $4-$byte) number. With the information provided, any image viewer that displays vector field data may be used to view the data. (See below for a different example where we have actually provided vectors in an Excel spreadsheet). The process by which these vectors were obtained is outlined below. Folder B consists of the same data now interpolated month by month resulting in a total of $51$ volumes over the time period of collection of the data. However the arrangement of the data in B is different. For the purposes of subsequent animation we arranged the data in a $256\times 256\times51$ volume which consists of a vector at each voxel in a single slice at the $51$ different times noted. As before, the vector is three 32-bit floating point numbers. Due to the size, we have separated these volumes into four sub-folders, comprising in sequence, the vectors at Slices $0-30, 31-70, 71-100,$ and $101-123$. Different views of the vectors as $51$ volumes in $256\times256\times124$ format are available, though they may also be reconstructed from the data in the repository.

We used the magnitudes of the vectors as images and constructed volume renderings (the opacity maps used for these renderings are available on request) of the $51$ such volumes. These are available in Folder C. With reference to the frame in which the data is constructed, the view point of the volume rendered maps were from a rotation around the $X-$axis of $30^{\circ}$, and around the $Y-$axis of $20^{\circ}$. Different viewpoints were selected according the viewing angle of rotation about the $Z-$axis as noted in the names of the files. Any particular frame of the animation, as well as the speed with which the animation is displayed, may be controlled through viewers available on the web. One such that we have used is available at https://onlineimagetools.com/gif-player as of February 15, 2021.

The process by which the vectors were obtained is now outlined. Some preprocessing of the data is required to obtain standard left-right alignment. A pullback transformation is applied to map a destination point to a source point (the first dataset). Thus there is a volume B which needs to matched to a source volume A, A and B being acquired adjacently in time, in that order. In our case both are registered to the coordinate system of the first data set which we call O. We do not describe the methods because such registration methods are described in detail in the source quoted below for the cortical surface data. Our procedure consisted of the following steps: (ⅰ) a voxel to scan coordinate transformation on A; (ⅱ) an inverse affine matrix to convert to the O space; (ⅲ) a local cross correlation optimized in each iteration step, and (ⅳ) minimizing a global mutual information measure. This allows us to generate a point correspondence between two datasets.

Three-dimensional vector fields on the cortical surface (Data sourced from ^[13] used for Figure 5). The cited paper may be consulted for the source and repositories of the original data at balsa@wustl.edu. As that website informs, BALSA stands for Brain Analysis Library of Spatial maps and Atlases and is a database for hosting and sharing neuroimaging and neuroanatomical datasets for human and primate species. The website also links to tools provided in the more general Connectome Project and Workbench for viewing and processing (via MatLab) the data. Using the Workbench, one can click on a "scene" file, which will open all files relevant to that figure and allows dataviewing. The surface files in Figure 3 scene are the ones we first used to create vector fields. The Excel file Garcia_ConsolidatedData_Large in the folder Garcia_et_al that we have provided is based on the surface files we mentioned, provided in GIFTI (.gii) format. The Connectome Workbench provides MatLab tools for reading such files, but we used our own Python programs to do so. In the different sheets of the file we have listed the surface points as well as the vector that leads from a point at one time to the corresponding one at the succeeding time.

However, we were informed by Dr. Kara Garcia that the surfaces referred to in the paper had been "resampled to the atlas resolution, which is higher than the resolution used for longitudinal registration and strain energy minimization." She provided new surface files "in which every vertex has been moved to the optimal position (rather than the higher resolution surfaces that contain lots of extra, unoptimized points).... Notably, each registration was run in both directions (older surface registered to younger surface, and vice versa), leading to two pairs of surfaces for each time period (folders ab, bc, cd, bd)." We therefore provide further Excel files noninjured.L.XY where XY = ab, bc, cd, and bd, respectively obtained by reading the coordinates from the surface files. The nomenclature of the surface files in each group used by Dr. Garcia and her collaborators is as follows:

"Forward" registration:

● YAS.LLR.surf.gii = original Younger Anatomical Surface mesh

● configincaltrelaxforward.anat.LLR.reg.surf.gii = older anatomical surface mesh after registration to younger

"Reverse" registration:

● OAS.LLR.surf.gii = original Older Anatomical Surface mesh

● configincaltrelaxreverse.anat.LLR.reg.surf.gii = younger anatomical surface mesh after registration to older."

Further we note the correspondence ab = Weeks $(27, 31)$; bc = Weeks $(31, 33)$; cd = Weeks $(33, 37)$ and bd = Weeks $(31, 37)$. We have provided the individual Excel files in self-explanatory nomenclature that were used to plot the vectors going from the earlier surface to the later. As stated, the plots in Figure 5 are vectors obtained from the Forward registration mentioned above.

Two-dimensional vector fields on a curve (Data referenced in ^[118], ^[119] and used in Figure 6). The data were provided by Dr. Shirley Bayer, originally in Adobe Illustrator form. The section images in the folder YakovlevDataAnnotatedPDF are the same but converted to pdf for easier viewing. Further we compiled them into volumes provided in the folder YakovlevData_VolumeSet which may be viewed again by any 3D Image viewing tool as for the first dataset above. All of the data in this set are $8-$bit (one byte). The data are taken from brain sections of fetuses at gestational weeks (GW) $29$, $31$, $34$, and $37$ which are not entirely dissimilar to the times when the above MRI datasets were obtained. The following information may be useful when viewing the volume data set. The data at GW $29$ is $550\times655\times52$, that at GW $31$ is $632\times 596\times91$, that at GW $34$ is $632\times596\times91$, and that at GW $37$ is $470\times582\times90$. The sections shown in Figure 6 are taken directly from the pdf's in the folder YakovlevDataAnnotatedPDF provided.

C.   Vector fields, Lie algebras and special functions

In this Appendix, we first mention some well known relationships between vector fields and Lie algebras, and on constructing Lie algebras from special functions. The latter material is taken from the extensive work of Feinsilver, Schott, and their colleagues, collectively hereafter referred to as FS. The specific papers consulted have been referenced in section 4, as is the book series by Feinsilver and Schott. Everything in these notes is to be found in their writings and we thank Professor Feinsilver for patiently explaining the work. We summarize very succintly their methods illustrated by the Heisenberg algebra. We hope that our description of the example shown is clear enough that it can be carried through without difficulty for any finite dimensional Lie algebra. We first summarize some introductory material on vector fields and Lie algebras.

C.1.   Vector fields and Lie algebras

The basic material on vector fields and Lie algebras that underlies the symmetry and equivalence approaches is widely available, and in partcular in ^[111]. A vector field is written in local coordinates as the first order differential operator or derivation

where ${\bf{x}} = \left(x^{1}, x^{2}, ....x^{d}\right)$ are local coordinates for the tangent space of a manifold. Such a vector field drives an evolution of points according to

where ${\bf{X}}$ is the $d-$tuple of the $X^{i}$'s, with some initial values specified. This is a one parameter flow, with parameter $t$. Given two vector fields in local coordinates in terms of functions ${\bf{X}}, {\bf{Y}}$, their Lie bracket is given by the commutator of the corresponding differential operators, and thus is also a vector field or, in local coordinates a derivation:

which is used to generate a Lie algebra by successive commutation. This is familiar from the "parallel parking" example. The generator $\frac{\partial}{\partial x}$ of a translation along $x$, the longitudinal axis of the car, and a rotation $y\frac{\partial }{\partial x}-x\frac{\partial}{\partial y}$about a vertical axis, allows one to translate along the transverse direction $y$. If the resulting Lie algebra is finite dimensional of dimensionality $N$, we have an $N-$parameter flow. Then, Lie's theorems when applicable allow one to state that smooth deformations of the manifold can be computed by composition of $N$ one-parameter flows. Mathematically precise statements of these remarks may be easily found, including in Olver's book ^[111]. Two different ways of writing an element of the group are used to develop the theory of generators of special functions. We therefore now turn to construction of these diffeomorphisms which are elements of the Lie group obtained by exponentiation of the Lie algebra. Only the structure of the Lie algebra is needed, not the specific representation in terms of the local coordinates, until we come to constructing special functions associated with representations on spaces of functions acting on the manifold in question.

C.2.   Transformations of coordinates

We now describe the FS constructions. However, we caution that we have made a few changes in notation so that the reader who cares to follow this while looking at the example in the papers can be warned. The changes are listed in Table 2. The assertions made without proof here are results that FS have obtained in various of their works, and are referenced in their papers. There are two fundamental constructions (from $X$ to $X^{L}$ and $X^{\prime}$, see below), related to left and right actions, and one fundamental result (from coordinates $\alpha$ of one kind (C.5) to coordinates $A$ of another (C.6) which is of the Baker-Campbell-Hausdorff type. Once these are understood, all the other assertions are easily stated and lead to various possibilities for construction of special functions. In this paper, and in most of their other work, FS use faithful matrix representations of the Lie algebras to carry through the first part of their computation to derive the ${X}^{L}$ and ${X}^{^{\prime}}$operators, as well as the associated $\pi$ matrices described below. However, since our work is confined here to representation spaces of functions of one real variable, we begin by using the conventional realization in terms of multiplication and differentiation operators. (Later, we indicate how to use matrix representations.) For reasons that will be clear we use $z$ to denote a real variable and argument of functions, as well as a multiplication operator by $z$ in the Lie algebra. No confusion should arise. We begin with the Heisenberg algebra in the form

with other commutators vanishing. We choose the representation $X_{1} =$ multiplication by $z$, $X_{2} =$ multiplication by a real number $h$ which then belongs to the center of the algebra, and $X_{3} = hd/dz$, all acting on functions of $z$. A group element may be represented either by

or by

where the $\alpha$ and the $a$ are scalars, and, without subscripts, $\alpha$ is meant to denote the sequence $\left(\alpha_{1}, \alpha_{2}, \alpha _{3}\right)$, and so for $a$. Without confusion, we also use $\alpha, a$ to stand for the corresponding list or vector in vector or matrix-vector equations, and sometimes for a single scalar. The context should make clear what is meant. FS call these coordinates of the first and second kind, respectively.

C.2.1.   Left and right vector fields

We define some differential operators. Write

The meaning of the arguments of $X_{i}^{{L}}$ will become clear in a moment. Arguments omitted in already defined functions should be understood to be present as when first introduced. We obtain these operators by pushing the $X_{i}$'s^* through $g$. Note that the expressions in (C.7) remain operators, implied to act on a function of $z$ to the right of the operator $g$. For example, we see that

^*We trust no confusion arises between our use of the apostrophe as here, and the primes which are used to denote various operators.

The other two operators cause no difficulty, so we have, regarded as a column vector,

The corresponding pi-matrix is defined by them as

with the notation

From these definitions, it follows that

Similarly, the following operator is defined, requiring pulling the Lie algebra operators through to the left:

Again the expressions are operators acting on a function of $z$ to the right of the expression in question. We have used a prime to denote a dual of the Lie algebra which is also a homomorphism thereof. Here, as is easily seen, the multiplication by $z$ is the only non-trivial operation:

The corresponding matrix satisfying

is thus seen to be

Defining the adjoint operator corresponding to an element $X$ on the Lie algebra vector space in the usual way as

we have

the last being true for our case since we have a $2$-step nilpotent algebra, so we may truncate the series at the second term. We then write

which yields the same results as equation (C.10). FS use faithful matrix representations to obtain the results. However, we need only the adjoint representation for this purpose, and this is very amenable to automating via a computer program. Recall that with the structure constants $c$ defined from

the $(kj)^{\rm{th }}$element of the adjoint matrix $X_{i} ^{\operatorname*{Ad}}$ representing $X_{i}$ is $c_{ij}^{k}$. Thus these matrices for (C.4) are

acting on the vectors corresponding to $X_{1}, X_{2}, X_{3}$

and which, since the algebra is nilpotent, cannot form a faithful representation. In this case, the adjoint action of both the zero of the Lie algebra and of $X_{2}$ are represented by the same matrix. However, the adjoint representation remains useful to obtain the operators $X^{L}, X^{\prime}$ since

where $\widehat{X}$ is the Lie algebra element corresponding to the vector argument. The definitions are meant to be perfectly general, not tied to this example, but it is clear how to proceed in more complicated cases. To summarize, we have

where the suffix $\pi_{i}^{\bullet}$ denotes the $i$'th row of the matrix. If $i = 1$ in equation (C.29) or $i = N$ in (C.30), we set the matrix in brackets to be the identity matrix $\mathcal{I}$. Alternatively, we can set the variable limit of the index to be $i$ in both cases instead of $i\pm1$ and perform the pointless operation $e^{a_{i} X_{i}^{\operatorname*{Ad}}}\chi_{i} = \chi_{i}$. In the Heisenberg algebra example, we see that

No summation of an infinite series is required to evaluate the exponential of any finite dimensional matrix representation of the algebra. However, the calculations soon become unwieldy enough that a symbolic manipulation program is most helpful. Further,

which follows immediately from the definitions. Note that the definitions of the $\pi-$matrices demand that the $X$'s in (C.32) be regarded as columns. We notice the following, which also happen to be general.

We now construct the double dual which will also be a homomorphism. Define, abstractly, $\mathcal{R}$, $\mathcal{V}$ (not usually adjoints of each other) to be raising and lowering operators for each element of the algebra satisfying:

and define $\mathcal{N}$ to be an ordering operator that moves all the raising operators to the leftmost position in an expression, and if there are lowering operators or derivatives, they are moved to the rightmost positions. Then

forms another homomorphism of the algebra. Explicitly for example, since the third row of $\pi^{L}\left(\mathcal{V}\right)$ would read $\left(0, \mathcal{V}_{1}, 1\right)$, we get

Thus

which satisfies the Lie algebra commutation.

The relation between the two kinds of coordinates is obtained by the method of characteristics, according to FS. Namely, consider the equations

where $\alpha$ is a row vector comprising the $\alpha_{i}$'s. Then

where $A\left(\alpha, 1\right)$ are the solutions at time $1$ of (C.39) with the integration constants set to zero. Henceforth, we denote $A\left(\alpha, 1\right)$ under these conditions simply as $A\left(\alpha\right)$. In this case we have,

Setting the integration constants to zero yields

We therefore have

The first set of canonical polynomials for the Lie algebra may be defined by considering $\overleftarrow{g}\left(\alpha, x\right)$ to be the generating functions of polynomials in the $x-$variables, where $x$ is a (set of) scalar(s). In other words, the canonical polynomials are $P_{n}\left(x\right)$ defined by the equality of the expressions below.

The righthand-most expression defines a multi-index notation for the integers and polynomials. Before we specialize to obtain some familiar polynomials, we introduce further quantities which may serve to define other canonical polynomial sets. First we define the Jacobian matrix

We have from (C.42)

Next we define

In other words we compute the Jacobian and replace the $\alpha$'$s$ by solving for them as a function of the $A$'s. Since $A\left(\alpha\right)$ is in general nonlinear in the $\alpha$, it is remarkable that this difficulty does not confront us in the Jacobian. Since $\alpha_{1} = A_{1}, \alpha_{3} = A_{3}$, we get, in terms of the values,

It is purely accidental that the functions $J, \overleftarrow{J}$ coincide in this case. Then a commuting set of differential operators in $x-$variables are defined from

where it is understood that the differential operators are to be placed rightmost in (C.49). As stated the $Y$'s are commuting operators, but they yield a remarkable result upon application to the constant $1$:

In other words $e^{\alpha_{1}Y_{1}+\alpha_{2}Y_{2}+\alpha_{3}Y_{3}}$acting on the constant $1$ is the function $\overleftarrow{g}\left(\alpha, X\right)$ given in (C.43). We shall assume $\mathcal{N-}$ordering of the differential operators in the expressions. In other words, we have

We can do a spot check on this. We find

where the $Y$'s are read off from (C.50). Speculatively, the joint eigenfunctions of Y may be worth exploring.

We now construct the third morphism. (FS denote both this and the second morphism above by the same symbols, but here we introduce an extra prime to distinguish this morphism from the previous ones). This is defined by

with

Computing the right hand side of (C.54),

and these operators form a morphism of our Lie algebra. Note that the relation between $\alpha$'s and $A$'s in

depends only on the commutation relations of the $X^{\prime}$s and so is valid for all the morphisms. For example, consider $X^{\prime\prime}$. We have from

Now let us apply both operators to the constant $1$^†. We can drop the lowering operator from the right hand side since it acts independently on the first variable. Then on the right hand side we are left only with commuting operators, so we may bring all the exponents together and add. We get

^†This may be considered the "vacuum" state from the point of view of quantum fields. However, note that the vacuum stated defined as what is annihilated by the " number" operator can be different. Thus $1$ is annihilated by $d/dx$, but the Heaviside function $H\left(x\right)$ is annihilated by $xd/dx$. $d/dx$ alone on $H$ will give the delta distribution which can be acted upon by $d/dx$ an infinite number of times, just as a constant can be acted upon by $x$ an infinite number of times.

We also note the following, which allows us to compute the matrix elements which are also special functions. From (C.49), let us posit

The action of raising and lowering operators on elements of the universal enveloping algebra are the familiar expressions from quantum theory, namely

Since the $Y$'$s$ commute, we can see from (C.51) that

where $e_{i}$ adds a $1$ to the power of the $i$'th coordinate function, and leaves all the others unchanged. Therefore, from the definition of the $x$'s in (C.63), we get the recursion relations

(see ^[121], below Eq (8) in their paper). Matrix elements of the group in this kind of number representation also form a class of special functions. Before we discuss this, we introduce yet another $\pi-$matrix, following ^[122]. We have already introduced notation in (C.55) for the transpose of the first $\pi -$matrix introduced, and also the matrix product that relates $X^{L}, X^{\prime}$, see (C.32). We therefore introduce

Then,

In words, the product of $\exp\left(aX_{i}^{\operatorname*{Ad}}\right)$, $a$ being a number and $X_{i}^{\operatorname*{Ad}}$ being the matrices in (C.26), result in the matrix $\bar{\pi}^{T}$. Our final formula from FS is one for the group elements. We use the multi-index notation, and use the Dirac notation $|n\rangle$ to represent an element

of the universal enveloping algebra. Then the matrix element $\mathcal{D} ^{nm}\left(a\right)$ is defined by

and the matrix elements are given as

the second line being the special case for an algebra consisting of three operators, in particular for the Heisenberg algebra we have been using as an example. Note that $X^{\prime}$ operators are expressed in terms of the $a$'s and their derivatives as indicated in (C.15) and worked out explicitly below that in equation (C.11). The $X^{\prime}$s on the right hand side of equation (C.74) are replaced by the multiplication and differentiation operators in the $a^{\prime}s$ and then applied to the monomials in the equation. The resulting expression is the relevant matrix element. These matrix elements satisfy recurrence relations when operated upon by the $X^{\prime}$, which result in further special functions, sometimes related to the hypergeometric functions, and are also "special functions" related to the algebra.

C.2.2.   Casimir operators

We have not made any use of these, but we quote the formula for the Casimir operators from ^[126] for a special case of Lie algebras, of the form illustrated in section 4 of the paper. Consider the algebra is generated by $q\left(D\right)$ and $x$ and has $n$ elements so that it is $n-1$–step nilpotent. Then the quadratic Casimir operators are of the form

We do not explore the eigenfunctions of such operators in part because they would not be useful for empirically derived vector fields. If $q\left(D\right) = D^{n}/n!$, these operators evaluate to zero. If we use the Lie algebra generators in the form $q\left(x\right), D$, the Casimir operators are multiplications by functions of $x$. ^[126] also exhibits formulas for Casimir operators higher than the quadratic.

C.2.3.   Summary

We may summarize the procedure as follows.

1. From the basic Lie algebra commutation table, construct the matrices of the adjoint representation.

2. Compute the $\pi^{L}$ and $\pi^{\prime}$ matrices from equations (C.29) and (C.30), respectively.

3. Obtain $A\left(\alpha\right)$ by solving equation (C.39).

4. Obtain the homomorphisms $X^{\prime}, X^{\prime\prime}, X^{\prime \prime\prime},$ and the commuting set $Y$ (see Table 2 for the defining equations).

5. Express the identity $f\left(\alpha, X\right) = g\left(A\left(\alpha\right), X\right)$ for the various morphisms.

6. From the previously expressed identities, search for special functions by acting upon suitable vectors such as the vacuum vector or eigenfunctions of the $Y$ operators.

7. Construct the matrix elements from the enveloping algebra and explore these for special functions.

C.3.   Special functions

We concisely list approaches to the construction of special functions suggested by FS. We did not make use of any of the methods advocated by FS in our paper. The particular nature of the Lie algebra allowed us to directly obtain the formulas displayed in section $4$ of the paper. Roads $3, 8,$ and $10$ were the ones illustrated in the paper.

C.3.1.   Ten Roads to special functions

We use $\alpha$ and $A\left(\alpha\right)$ to mean the respective lists, and we specialize to the case where the generators of the Lie algebra are written $q\left(x\right), D$, though the methods can be generalized to other Lie algebras. Further, $\alpha_{i}$ is associated with $D^{i-1} q(x), i = 1, 2, ...$ and in this note we will associate $\alpha_{0}$ with $D$. $X$ will be any morphism of the Lie algebra, and the indexing will be the same, i.e., $X_{0}\sim D, \; X_{1}\sim q$ and so on. The order will be $A_{1}, A_{2}, .....A_{N}, A_{0}$ where the Lie algebra is $N-$step nilpotent.

Road 1: Generating functions Apply ${ \prod} \exp\left(A_{i}X_{i}\right)$ to some function $f\left(x\right)$, conventionally chosen to be $1$. Then consider the polynomials

where $J$ is a multi-index.

Road 2. Generating functions, commuting operators We can apply the same method as above but this time where the $X^{\prime}$s are replaced by the commuting differential operators $Y$ in the above.

Road 3. Orthogonal polynomials $q(x)-D$ and $D$ are operators that are formal adjoints of each other in the space of functions which are square integrable with respect to the weight function $\exp\left(-Q\left(x\right) \right)$, where $Q^{\prime}\left(x\right) = q\left(x\right)$. They can be made adjoints with suitable zero boundary conditions for the functions, or with addition of delta functions to address non-zero boundary conditions. It is then natural to consider polynomials orthogonal with respect to the inner product suitably defined with this weight. Let $p_{n}\left(x\right)$ be such polynomials.

Road 4. Fourier transforms The Fourier transforms

where $\mathcal{F}$ denotes the Fourier transform, can be checked to form an orthogonal family of functions of $k$ as well.

Road 5. Canonical transform Since $q\left(-D\right), x$ generate the same Lie algebra as $q\left(x\right), D$, we may evaluate ${ \prod} \exp\left(A_{i}X_{i}\right)$ without the need to apply it to any function and follow the routes given in $1$ and $2$ above.

Road 6. Specialization of the $\alpha^{\prime}$ s In the above differentiations, $J$ was a multi-index since there were many $\alpha^{\prime }$s. We may reduce the number of these in various ways. One way is to scale $\alpha_{0} = \lambda = \alpha_{1}$ and the other $\alpha^{\prime}$s scale as $\lambda^{i}$ so that each $A_{i}$ becomes a homogeneous polynomial in $\lambda$. Then the product of exponentials is a function of the single parameter $\lambda$ and we may differentiate with respect to this. This seems the natural choice, but another is to choose the scale for $\alpha_{0}$ to be a constant (i.e., $\alpha_{0} = \lambda^{0}$) in which case we can choose all the other $\alpha^{\prime}$s to scale linearly with $\lambda$.

Road 7. Matrix elements Illustrated in the discussion following equation (C.74).

Road 8. Eigenfunctions I We may square the Lie algebra elements and seek the eigenfunctions of $q(x)^{2}+D^{2}$ or similar self adjoint operators in a suitable Hilbert space of functions.

Road 9. Eigenfunctions II From the remarks in Road 3, $\left(q\left(x\right) -D\right) D$ is also self adjoint in a suitable Hilbert space. Find the eigenfunctions of these.

Road 10. Application of generator Apply $q\left(x\right) -D$ successively to a suitable function (the constant function $1$ will do), and orthonormalize the resulting polynomials or functions. We have discussed this in section 4, and we make the following additional remarks. We see therefore that

can be taken as a generating function for special polynomials $P_{n}$. In the case $Q = x^{2}$, we get the generating function for the Hermite polynomials. In general by differentiating with respect to $\alpha$, we get

By expanding $q$ in powers of $\alpha$, and equating coefficients of like powers of $\alpha$, we get a recursion relation in the index of the polynomials which has as many terms as there are non-zero differential coefficients of $Q$. So for the quadratic, we recover the familiar three term recurrence for the Hermite polynomials. Similarly if we differentiate with respect to $x$, we get

Again, expansion of $q\left(x-\alpha\right)$ in powers of $\alpha$ will yield a recurrence expressing $P_{n}$ in a series with terms which are derivatives of lower order, and by inversion of the series, $P_{n}^{\prime}$ may be expressed in terms of the polynomials of lower order. These expressions generalize recursion relations obtained for the Hermite polynomials, but orthogonality is lost. However, any finite set of polynomials may be orthogonalized via the Gram-Schmidt or other procedure.