Quasicrystal patterns in a neural field model

Aytül Gökçe ,1,* Stephen Coombes ,2,† and Daniele Avitabile 3,4,‡ 1Department of Mathematics, Faculty of Arts and Science, Ordu University, 52200 Ordu, Turkey 2Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, NG7 2RD, United Kingdom 3Department of Mathematics, Vrije Universiteit Amsterdam, Faculteit der Exacte Wetenschappen, De Boelelaan 1081a, 1081 HV Amsterdam, Netherlands 4Mathneuro Team, Inria Sophia Antipolis, 2004 Rue des Lucioles, 06902 Sophia Antipolis Cedex, France

(Received 12 December 2019; accepted 3 February 2020; published 2 March 2020) Doubly periodic patterns in planar neural field models have been extensively studied since the 1970s for their role in explaining geometric visual hallucinations. The study of activity patterns that lack translation invariance has received little, if any, attention. Here we show that a scalar neural field model with a translationally invariant kernel can support quasicrystal solutions and that these can be understood using many of the theoretical tools developed previously for materials science. Our approach is constructive in that we consider constraints on the nonlocal kernel describing interactions in the neural field that lead to the simultaneous excitation of two periodic spatial patterns with incommensurate wavelengths. The resulting kernel has a shape that is a modulation of a Mexican-hat kernel. In the neighborhood of the degenerate bifurcation of a homogeneous steady state, we use a Fourier amplitude approach to determine the value of a Lyapunov functional for various periodic and quasicrystal states. For some values of the parameters defining a translationally invariant synaptic kernel of the model, we find that quasicrystal states have the lowest value of the Lyapunov functional. We observe patterns of 12-fold, 10-fold, and 6-fold rotational symmetry that are stable, but none with 8-fold symmetry. We describe some of the visual hallucination patterns that would be perceived from these quasicrystal cortical patterns, making use of the well known inverse retinocortical map from visual neuroscience.

I. INTRODUCTION
One of the most well known success stories in mathematical neuroscience is the theory of geometric visual hallucinations, developed by Ermentrout and Cowan [1] and subsequently extended by Bressloff et al. [2]. These authors viewed the cortex as a planar pattern-forming system that could be described by a neural field equation. Neural fields are typically cast in the form of integro-differential equations with a nonlocal spatial interaction term that is specified by a synaptic kernel function describing anatomical connectivity. Although at first sight the spatially nonlocal evolution equation for a neural field is somewhat different from local partial differential equation (PDE) models of pattern-forming systems, such as the Swift-Hohenberg equation, it can be analyzed with many of the same tools. These include linear Turing instability analysis, weakly nonlinear analysis, and symmetric bifurcation theory to name just the most common ones. These have been put to use to determine what geometric visual hallucinations can tell us about the primary visual cortex (V1) [3]. The key observation of Ermentrout and Cowan was that the geometric shapes, originally classified by Klüver, in terms of form constants for (i) tunnels and funnels, (ii) spirals, (iii) lattices including honeycombs and triangles, and (iv) cobwebs, could all be described in terms of periodic patterns of activity occurring in V1. Moreover, they attributed the emergence of such patterns to a Turing instability of a homogeneous steady state, arising through the variation of a control parameter representing the effect of a hallucinogenic drug such as LSD, cannabis, or mescaline on the cortex. Thus the problem of understanding hallucinations was reduced to finding the kernel of a neural field model that would support the relevant type of Turing instability. Beyond bifurcation the emergence of striped, square, or hexagonal patterns of neural activity in V1 induce visual percepts, which can be understood with further information about how retinal activity is mapped to cortical activity. As a first approximation (away from the fovea) this map is conformal and takes the form of a complex logarithm [4]. When applied to oblique stripes of neural activity in V1 the inverse retinocortical map generates hallucinations comprising spirals, circles, and rays, while lattices like honeycombs or checkerboards correspond to hexagonal activity patterns in V1. Ermentrout and Cowan found that patterns of translationally invariant synaptic interaction with short-range excitation and long-range lateral inhibition were the answer. In reduced versions of their neural field model that consider only a scalar field (rather than a pair of coupled fields) on an infinite plane, the effective interaction kernel can be described using a Mexican-hat function. Somewhat surprisingly, the search for kernels that can support nonperiodic patterns has received very little, if any, attention. We say surprisingly because of the interest in the condensed matter communities of finding PDE models that support such solutions [5]. This has mainly been motivated by the discovery of physical quasicrystals [6], in which atoms are ordered over long distances but not in the periodically repeating arrangement of traditional crystals, for which the Nobel Prize in chemistry was awarded to Shechtman in 2011. Dodecagonal quasicrystals in soft condensed matter systems have also been found experimentally [7,8]. For a very recent discussion of how multiple spatial scales can lead to the formation of quasicrystals in soft condensed matter systems, we refer the reader to [9].
Of course, a great deal of interest in quasicrystals has also been motivated by mathematics, specifically as it relates to tiling, with perhaps the most well known result being due to Penrose for his work in the 1970s on a set of just two tiles that can tile the plane. These aperiodic tilings were generalized to three dimensions by Levine and Steinhardt [10] and used to discuss the structural properties of three-dimensional quasicrystals, including their diffraction patterns [10,11]. From a pattern-forming perspective, it has been recognised for some time that one way to make models that support quasicrystal solutions is to consider nonlinearly interacting systems with two incommensurate spatial scales, as in the Landau model of Mermin and Troian for quasicrystalline order [5] or the three-wave (nonlinear) interaction exemplified by the phenomenological PDE model of Rucklidge et al. [12] for describing quasicrystal patterns seen in the Faraday wave experiment. In these experiments strong vertical vibrations of a viscous liquid layer can elicit spatially complex standingwave quasipatterns at parameters where two length scales in the correct ratio are excited or weakly damped; see work by Skeldon and Rucklidge for a discussion of their stability [13]. In this paper we build on these general principles to show how to construct a neural field model that can support stable quasicrystal patterns of activity.
In Sec. II we describe the simple scalar two-dimensional neural field model that is the subject of this paper. We also present a Lyapunov functional for the model. The nonlocality of the model is encoded in an anatomical interaction kernel and is the subject of Sec. III. Here we introduce a set of conditions for a translationally invariant kernel to support a Turing instability that will excite two incommensurate spatial scales, whose nonlinear interaction may lead to a quasicrystal pattern. This nonlinearity is represented in the model with a two-parameter sigmoidal firing rate function. In Sec. IV we use the Lyapunov functional to construct a phase diagram for the model showing the most stable states corresponding to dodecagonal (DDQC), decagonal (DQC), and octagonal (OQC) quasicrystal solutions, as well as hexagonal (HEX) periodic patterns. For some choices of the sigmoidal firing rate and synaptic kernel parameters we find that quasicrystal solutions have the lowest value of the Lyapunov functional. This is confirmed with direct numerical simulations of the neural field model. The perception of these aperiodic patterns is described in Sec. V, making use of the known inverse retinocortical map from visual neuroscience. Finally, in Sec. VI we discuss the shape of the kernel that gives rise to quasicrystal neural activity patterns and suggest that it may be realized in V1 via a spatially varying modulation of anatomical connection strengths.

II. NEURAL FIELD MODEL
Neural field models are continuum descriptions of neural tissue, which aim to capture coarse-grained dynamics for synaptic neural activity. Since their inception in the 1970s, they have been used in a variety of neural contexts, to explain phenomena ranging from visual hallucinations to brain rhythms. For a recent collection of survey articles about their use we refer the reader to [14]. Here we will consider one of the simplest variants of a neural field, namely, one that describes a single neural population over a flat planar domain. Introducing the scalar quantity u = u(r, t ), with r ∈ R 2 and t 0, the model takes the form of an integro-differential equation, which we write in the compact form Here the symbol ⊗ denotes spatial convolution such that The anatomical kernel w is given by a function of Euclidean distance and the nonlinearity f represents a transduction rule for converting the neural activity u into a firing rate.

Lyapunov energy functional
For the case that the firing rate f is monotonically increasing and the kernel w is even and decays sufficiently quickly away from the origin, there is a Lyapunov functional for (1) [15]. We will assume these properties from now on. In particular, we will take the firing rate to be a sigmoidal function that depends on two parameters h and β. The former is interpreted as a threshold value (describing when f switches from a low to a high value) and the latter controls the steepness of the sigmoid at threshold. We note that for this choice we have the Riccati relation Figure 1 shows results of a direct numerical simulation of an evolving decagonal quasicrystal with a specific choice of synaptic kernel w. This example nicely illustrates the ability of neural field models to generate quasiperiodic patterns from spatially localized initial data. For further details see the discussion around Eq. (9), as well as Sec. III (for the synaptic kernel) and Appendix E for details of the numerical scheme. The Lyapunov functional can be written in the form It is straightforward to establish that the Lyapunov functional is decreasing along trajectories such that The mean-value theorem can also be used to calculate the functional derivative of E with respect to u as (6) for some smooth test function φ = φ(r, t ). Thus the critical points of E [u] satisfy (−w ⊗ f (u) + u) f (u) = 0, and since f (u) = 0 this recovers the steady state of the model (1). In Fig. 2 we plot the Lyapunov functional from Eq. (4) throughout the evolution of the growing decagonal quasicrystal, illustrated in Fig. 1. Each layer of growth around a central localized quasipattern is accompanied by a gradual decrease in the Lyapunov functional. This also shows that a global decagonal quasicrystal (covering the whole plane) tends to a phase with the lowest value of the Lyapunov functional. In a generic situation we have that the Lyapunov functional of a solution trajectory decreases as t → ∞, critical points of the Lyapunov functional are steady-state solutions, and minima are stable. Thus the Lyapunov functional (4) is well suited to determining the stability of both spatially periodic and quasicrystal solutions. The former are well known to exist in neural field models with kernels exhibiting short-range excitation and long-range inhibition [1]. The conditions for the latter to exist in a neural field model do not appear to have been previously explored, although it is well known for many physical systems that nonlinear interactions between patterns at two length scales can stabilize quasicrystalline order. We now develop a simple Turing analysis to determine the type of kernel function that can give rise to a pattern with two excited wavelengths.

III. KERNEL WITH TWO SPATIAL SCALES
In a neuroscience context it is very common to meet tissue models that respect a balance between excitation and inhibition. This is readily imposed on the model given by Eq. (1) by enforcing the condition that R 2 w(|r|)dr ≡ W = 0. Although not strictly necessary for the discussion that follows, it does simplify some of the analysis and in particular enforces a unique fixed point. To see this, consider a space-and timeindependent homogeneous steady state given by u(r, t ) = u 0 for all r and t. Substitution into (1) gives a nonlinear equation for u 0 as u 0 = f (u 0 )W . Thus, for the choice W = 0 we can sidestep the need to find the roots of a nonlinear equation and have simply that u 0 = 0. The evolution of perturbations around u 0 can be determined by substitution of u(r, t ) = u 0 + v(r, t ), for some small v, into Eq. (1) and expanding to first order. The linearized equation describing the evolution of perturbations has a space-time separable solution v(r, t ) ∼ e λt e ik·r , with where γ = f (0) andŵ(k) is the Fourier transform of w, Since the kernel is an even function w(r) = w(|r|), then w(k) =ŵ(|k|), and λ given by Eq. (7) is a real function of |k|. If the Fourier transform of w has a maxima away from the origin, say, at k = |k| = k c > 0, then a Turing instability to a spatially periodic pattern, with wavelength 2π/k c , will occur when 1/γ decreases throughŵ max =ŵ(k c ). Moreover, ifŵ max is doubly degenerate in the sense thatŵ(1) =ŵ max = w(q) for two distinct wave numbers k = 1 and k = q, then the Turing instability analysis would predict the excitation of a pattern with two distinct spatial scales in the neighborhood of the bifurcation. This suggests a constructive approach to find a kernel that would support the generation of quasicrystal patterns in the model (1), namely, to consider the shape of the Fourier transform of the kernel function w to enforce both balance and degeneracy. The former condition is described bŷ w(0) = 0 and the latter byŵ(1) =ŵ max =ŵ(q),ŵ (k)| k=1 = 0 =ŵ (k)| k=q , andŵ (k)| k=1,q < 0. The parameter q is an irrational number that ensures that the two excited wave numbers are incommensurate. As a candidate choice for w let us consider a kernel model for short-range excitation and long-range inhibition, but with some degree of modulation. For ease of computing Fourier transforms we make the explicit radially symmetric choice where After introducing the function we may writeŵ(k) as In Appendix A we show that A numerical optimization can now be used to find the parameter values of α 1 , α 2 , b 1 , b 2 , s 1 , and s 2 to obtain the required shape ofŵ(k). Here the parameters b 1 and b 2 are set in such a way thatŵ(1) andŵ(q) are numerically extremely close. Then five simultaneous equationsŵ (k)| k=1,q = 0,ŵ(0) = 0, andŵ(k)| k=1,q −ŵ max = 0 are solved for the parameters α 1 , α 2 , s 1 , s 2 , andŵ max so that the balance and degeneracy conditions for the kernel (9) can be satisfied for particular values of α 1 , α 2 , b 1 , b 2 , s 1 , and s 2 . An example is shown in Fig. 3 for the choice q = 2 cos(π/5), demonstrating the shape of the kernel in real and Fourier space. The function w plotted in Fig. 3(a) describes how neurons a distance r apart influence one another and is in the form of a Mexican hat with a modulated tail. In Fig. 3(b) we also show a zoom around the two wave numbers k = 1 and k = q illustrating the degeneracy in the maxima ofŵ(k) such that a quasicrystal with two length scales can be excited whenŵ max = 1/γ . In Fig. 4  It should be noted, from the power spectra, that the two basic length scales of the model are determined from wave vectors lying in the annuli centered around the circles k = 1 and k = q. The type of emergent pattern is determined with the proper choice of the second length scale q, for which q = 2 cos(π/12), q = 2 cos(π/5), and q = 2 cos(π/8) are considered for dodecagonal [ Fig. 4(a)], decagonal [ Fig. 4(b)], and octagonal [ Fig. 4(c)] quasicrystal patterns, respectively. The real space and Fourier space of these structures present 12-fold, 10-fold, and 8-fold rotational symmetries. Here dodecagonal and decagonal quasicrystal patterns evolve from a localized perturbed spot, while the octagonal quasicrystal is initiated with an 8-fold rotationally symmetric pattern. Initial data for each candidate structure are given in Appendix D by Eq. (D1). Note that the initial configuration for an octagonal structure with a desired 8-fold symmetry ultimately converges to other states, more commonly hexagonal, though its Lyapunov functional evolves slowly. Figure 4 including dodecagonal, decagonal, and hexagonal patterns [16,17].

IV. FOURIER AMPLITUDE DESCRIPTION
The stability of quasicrystals is expected to be promoted by nonlinear wave interactions between three or more waves (see [12] for further discussion). In a physical context, for either hard materials such as metallic alloys or soft condensed matter systems, the use of a free-energy functional has proven especially useful in clarifying the role of nonlinear interactions. A case in point is the free energy of the Lifshitz-Petrich model [16,18], itself a modified version of the Swift-Hohenberg model constructed to develop quasicrystalline patterns, for which decagonal and dodecagonal quasicrystals can be stable whereas the octagonal quasicrystal state is metastable. Thus, given that the neural field model (1) possesses a Lyapunov functional given by (4), we may naturally employ techniques from soft condensed matter systems to study phase diagrams of the system. In particular, the Lyapunov functional (4) has much in common with that used to describe block copolymers [19], for which a cut-and-project method can be used to calculate the free-energy functional of quasicrystals. This method was originally proposed by Meyer [20] and projects the higher-dimensional periodic lattice points within a stripe onto a lower-dimensional space to obtain a quasicrystal. For a recent discussion of this approach in the context of soft condensed matter physics we refer the reader to the recent work of Jiang et al. [17]. They show that with this approach the computation of the free energy can be readily carried out in Fourier space and gives a unified framework for studying both periodic crystals and quasicrystals. To exploit the cut-and-project methodology we first consider a useful approximation of our Lyapunov functional valid in the neighborhood of an instability of the homogeneous steady state.

A. Logit approximation
The sigmoidal firing rate function has a logit function inverse g = f −1 , where In the neighborhood of u = a we may expand g(u) as a Taylor series where η 0 = g(a), η 1 = g (a), η 2 = g (a)/2, and η 3 = g (a)/6. A simple calculation gives g It is now convenient to introduce z(r, t ) = f (u(r, t )) and write (4) as Under a change of variables the second term on the right-hand side of (15) can be written as so that if z(r, t) remains in the neighborhood of f (0) we may use the expansion of the logit function to obtain where we use the result that g( f (0)) = 0 (so that in this case η 0 = 0). For a balanced kernel we have that E [z 0 ] = 0, with a homogeneous steady state z 0 = f (0).

B. Projection method
From a mathematical perspective, quasicrystals can be viewed as projections of a higher-dimensional lattice. To be more specific, the Fourier transform of a d-dimensional quasicrystal is nonzero only at a dense set of points spanned by integer multiples of a finite set of basis vectors (which are the projections of the primitive reciprocal lattice vectors) of an n-dimensional lattice, with n d [21]. The n-dimensional reciprocal vectors can be spanned by a set of bases b i , which are the primitive reciprocal vectors in the n-dimensional reciprocal space, with integer coefficients. Hence, an n-dimensional reciprocal vector H can be written as H = n i=1 h i b i ∈ R n , with h i ∈ Z and b i ∈ R n (canonical basis vectors). Introducing a d × n projection matrix S allows the construction of a physical d-dimensional wave vector k as k = SH. Here n and S are chosen in a problem-specific manner, and projections for generating 12-fold, 10-fold, and 8-fold symmetric quasicrystals are given in Appendix B. Using the n-dimensional periodic lattice and the projection matrix S, the expansion of any d-dimensional quasiperiodic function ψ (r), r ∈ R d , can be written where B = (b 1 , . . . , b n ) ∈ R n×n and h = (h 1 , . . . , h n ) ∈ Z n . The terms φ(h) are recognized as Fourier coefficients, though it is important to remember that (18) is a d-dimensional quasiperiodic function determined as a slice of an ndimensional periodic structure whose orientation is determined by S. It is natural to assume that close to a bifurcation from a homogeneous steady we may write where k is a vector generated from a high-dimensional periodic lattice according to k = n i=1 h i Sb i and 1. Using (19), we may now evaluate the Lyapunov functional given by (17) close to a bifurcation, in terms of the Fourier amplitudes φ(h). Choosing a scaling of parameters such that γ −1 ∼ O( 2 ) ∼ η 1 and η 2 ∼ O( ), then all terms in (17) contribute equally. Considering a set of new parameters E = 4Ẽ , w = 2ŵ , η 1 = 2η 1 , η 2 = η 2 , and η 3 =η 3 , the Lyapunov With a minor modification of the approach used by Jiang et al. [16], we may now easily compute the Lyapunov functional for various quasicrystal and periodic patterns in terms of the Fourier amplitudes on the circles of |h| = 1 and |h| = q, respectively. Denoting these amplitudes by φ 1 and φ q , respectively, we may then minimizeẼ =Ẽ (φ 1 , φ q ) for a given pattern type. The forms of the Lyapunov functional for dodecagonal, decagonal, and octagonal quasicrystals are given in Appendix C, together with that of a periodic hexagonal pattern. By comparing the minimum of the Lyapunov functional for each pattern configuration at a given point in parameter space we may also determine a phase diagram. Here we do this for the two scaled parameters (η 1 ,η 2 ) (see Appendix C). For fixed values of α 1 , α 2 , b 1 , b 2 , s 1 , s 2 , and q in the synaptic kernel, the phase diagram over a range of η 1,2 values is shown in Fig. 5, where the most stable states for DQC and HEX patterns are shown as an example. In the black region, e.g.,η 1 > 50.63, all patterns mentioned above have the same (zero) energy [using the approximation given by (17)], as does the homogeneous steady state. In simulations of the full nonlinear model we typically see the evolution to the homogenous steady state in the black region. Asη 1 is decreased, the decagonal structure is favored, and with a further decrease, e.g.,η 1 < 46.9, the hexagonal phase has the FIG. 6. Perception of (a) dodecagonal, (b) decagonal, and (c) octagonal quasicrystals patterns of neural activity in V1, as presented in Fig. 4, constructed via the inverse retinocortical map. lowest value of the Lyapunov functional for allη 2 values. It can be seen from Fig. 5 that the decagonal pattern is stable over a narrow range ofη 1 for large values ofη 2 . Considering the parameter setting of the kernel for 10-fold quasicrystals in Fig. 4(b), stable dodecagonal and octagonal quasicrystal patterns are not found for any choice ofη 1 andη 2 values in Fig. 5. Repeating the energy comparison among the candidate structures with parameters from the 12-fold pattern in Fig. 4(a), only hexagonal and dodecagonal symmetric patterns are observed (not shown). A stable octagonal pattern is not found for any choice of parameters. Note that an increase in theη 2 value corresponds to a decrease in h and β values in the sigmoidal firing rate function (3), whereas larger values ofη 1 correspond to smaller β and larger h.

V. VISUAL HALLUCINATIONS
One of the main structures of the visual cortex is that of retinotopy, a neurophysiological projection of the retina to the visual cortex. The log-polar mapping [4] is perhaps the most common representation of the mapping of points from the retina to the visual cortex and is key to understanding some of the visual hallucinations that can be induced by hallucinogenic drugs and in particular their geometry. These can be categorized by the Klüver form constants: tunnels and funnels, spirals, lattices including honeycombs and triangles, and cobwebs. It was the great insight of Ermentrout and Cowan [1] that these could be recovered after an application of the inverse retinocortical map to spatially periodic activity arising from a Turing instability in V1. The action of the retinocortical map turns a circle of radius r in the visual field into a vertical stripe at x = ln(r) in the cortex and also turns a ray emanating from the origin with an angle θ into a horizontal stripe at y = θ . Simply put, if a point on the visual field is described by (r, θ ) in polar coordinates, the corresponding point in V1 has Cartesian coordinates (x, y) = (ln(r), θ ). Thus to answer how a quasicrystal pattern would be perceived we need only apply the inverse log-polar mapping. Figure 6 shows how neural activity patterns for dodecagonal, decagonal, and octagonal quasicrystals in V1 are mapped to various exotic percepts in the visual field and extend slightly the class of visual hallucinations that can be captured by neural field models.

VI. DISCUSSION
In this paper we have considered a simple scalar twodimensional neural field model and shown how one can construct a translationally invariant anatomical interaction kernel to guarantee the emergence of quasicrystalline activity patterns beyond a Turing instability. In line with current thinking of how quasipatterns can be generated in physical systems, this requires an interaction between two incommensurate spatial scales whose ratio is given by the parameter q in the synaptic kernel. We have shown that one such kernel that achieves this is a form of modulated Mexican hat. To determine the stability of emergent quasiperiodic patterns we have made use of the existence of a Lyapunov functional for the neural field model. This has allowed us to tap into approaches previously deployed for theoretical studies of condensed matter systems exhibiting quasicrystal structure and in particular those for the determination of phase diagrams (via minimizing the Lyapunov functional of the system). Similarly to many soft condensed matter systems, we have found there are regions in parameter space that support stable dodecagonal and decagonal quasicrystals, while octagonal quasicrystals invariably have a higher value of the Lyapunov functional and are metastable.
The modulated Mexican hat has been shown to robustly support stable quasicrystals with an n-fold rotational symmetry with n = 10, 12 (and an 8-fold solution only being metastable). However, in principle, there is nothing to preclude the construction of quasipatterns with other degrees of rotational symmetry. The projection method that we have used here, for so-called rank-4 quasicrystals (describing 8-fold, 10-fold, and 12-fold symmetric quasicrystals projected from a four-dimensional space), can also be deployed to generate other n-fold symmetric quasicrystals by projection from higher-dimensional spaces, for example, 7-fold, 9-fold, and 18-fold symmetric quasicrystals from a six-dimensional space. However, it is worth noting that in many physical systems there are often geometrical constraints which impede the formation of quasicrystals with other than 5-fold or 10fold symmetry (see, e.g., [22]).
Interestingly, kernels described by periodically modulated exponentials have previously been invoked as models for interaction in the prefrontal cortex [23] (motivated by labeling studies showing that approximate periodic stripes are formed by coupled groups of neurons). Moreover, given the known periodic modulation of primary visual cortex on the millimeter scale, as reviewed in some detail by Bressloff [24], it is not unreasonable that a modulated Mexican-hat shape is more representative of interactions in primary visual cortex than an unmodulated one. The question about whether quasicrystal patterns are in fact ever observed in drug-induced hallucinations or as a result of certain classes of visual stimuli is a fascinating one, worthy of further attention by the psychophysics community. The work presented here has focused on the spontaneous generation of patterns and not ones induced directly by sensory stimuli. However, given that quasipatterns in physical (Faraday) systems can be excited by periodic temporal forcing [25], this motivates a protocol for human studies. Indeed, flicker-induced hallucinations were previously studied from a theoretical perspective in neural fields with time periodic forcing by Rule et al. [26], and it would be very natural to extend the work here to include models of spatiotemporal sensory drive.
Looking forward, it would be interesting to use techniques described by Rankin et al. [27], for numerical bifurcation and continuation for homoclinic snaking, to better understand how quasipatterns can grow by adding structure to the outside of a central localized quasipattern, as illustrated in Fig. 1. Finally, it is worth mentioning that the techniques presented here could also be developed to study superlattice patterns that display periodic spatial structures having two separate scales [28].

APPENDIX A
Here we show how to evaluate the integral representation of the function F (k, a, σ ) given by (11). Using plane polar coordinates, we have that The last integral can be evaluated with the change of variable z = e iθ , performing a contour integral around the unit circle, and noting that there is a single contribution from a pole at i(m − √ k 2 + m 2 )/k. This gives the result (12).

APPENDIX C
Here φ 1 and φ q are real amplitudes on the circles of |h| = 1 and |h| = q, respectively. We use the result that vectors on the |h| = q circle are sums of two neighboring reciprocal vectors on the |h| = 1 circle. Absorbing a factor of (2π ) 2 within the energy functional, we adapt the calculations in [16] (with minor modification) to find where φ 1 and φ q are real amplitudes on the circles of |h| = 1 and |h| = q, for which q = 2 cos(π/12), q = 2 cos(π/5), and q = 2 cos(π/8) respectively determine the type of pattern.
Here the label HEX refers to a hexagonal periodic pattern.

APPENDIX D
Initial data for the neural field model (1) are chosen in the form u(r, 0) = e −a|r| 2 /L N j=1 (e ik j ·r + e iq j ·r ), k j = cos θ j sin θ j .
(D1) Here the vectors k j , j = 1, . . . , N, are equally spaced around the circle with |k| = 1 and we set q j = qk j . Figure 1(a) shows the initial condition for a 10-fold quasicrystal with a = 2, L = 36π , N = 10, and q = 2 cos(π/5), by which a localized spot is perturbed with a 10-fold quasiperiodic pattern. Figures 1(b)-1(e) show a transient state after 40, 100, 120, and 160 time units, respectively, and the equilibrium state after 300 time units. Note that we set a = 0 for an octagonal quasicrystal so that the domain is initially filled with an 8fold pattern with L = 30π . To visualize quasicrystal patterns better, the domains in Figs. 1, 4, and 6 have been truncated to [−24π, 24π ] × [−24π, 24π ].

APPENDIX E
The numerical simulation of the full model (1) was performed in the plane by discretizing (1) in space on a regular square mesh and solving the resultant set of ordinary differential equations using MATLAB. In our simulations, a pseudospectral evaluation of the convolution w ⊗ f (u) is performed using a fast Fourier transform (FFT), followed by an inverse FFT on a large square computational domain. The Fourier transform of w ⊗ f (u) takes the product form w(|k|) f (k), and this provides substantial computational speedup over quadrature-based numerical methods for calculating w ⊗ f (u). We set a grid of N = 2 10 equally spaced points along each spatial dimension and used MATLAB's inbuilt ODE45 algorithm to evolve the system forward in time.
Quasicrystals are ordered structures that fill the space without spatial periodicity. Hence, the main limitation of using the FFT algorithm is that small errors may arise from the choice of a periodic computational domain. To overcome this limitation and enable the efficient use of the FFT algorithm, direct numerical simulations are performed on a large domain.