Spatiotemporal Kernel of a Three-Component Differential Equation Model with Self-control Mechanism in Vision

This paper examines a three-component differential equation model with a self-control mechanism in vision as a slight extension of the lateral inhibition model proposed by Peskin (Partial differential equations in biology: Courant Institute of Mathematical Sciences Lecture Notes, New York, 1976). First, we derive a condition under which the exact solution of our differential equation model for time-dependent input I=I(t)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I=I(t)$$\end{document} is described by the convolution integral with a temporal biphasic function. Second, we analyze the model with the input signal I=I(t,x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I=I(t,x)$$\end{document} depending on time t and position x∈R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x \in \textbf{R}$$\end{document}, and we prove that the solution can be represented in convolution integral form and that t1>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_1>0$$\end{document} exists such that the spatiotemporal integral kernel Ku(t1,x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u(t_1,x)$$\end{document} is positive for x∈R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x \in \textbf{R}$$\end{document} and t∈(0,t1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t \in (0,t_1)$$\end{document}. Moreover, we numerically demonstrate that there exists t2(>t1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_2 (> t_1)$$\end{document} such that Ku(t2,x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u(t_2,x)$$\end{document} includes the Mexican-hat function and a temporal biphasic function under certain parameter conditions. From these mathematical and numerical results, we find that there is a time lag before the Mexican-hat function appears in Ku(t,x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u(t,x)$$\end{document}, and the shape of Ku(t,x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u(t,x)$$\end{document} is similar to the receptive field structure observed in experiments in the field of neurophysiology. We conclude that the partial differential equations for visual information processing can be used to analytically determine the shape of the spatiotemporal kernel indicating the self-control mechanism.


Introduction
When we observe an object, visual information processing occurs not only in the brain but also in the retina.It is well known that the retina plays an important role in visual information processing, including color vision and adaptation [4,48,52].The retina is a hierarchical cellular tissue constructed by principal nerve cells, such as photoreceptor cells, horizontal cells, and retinal ganglion cells [16].The photoreceptor cells convert input stimuli to electrical signals, and the electrical signals are transmitted to the retinal ganglion cells by receiving feedback from the other nerve cells.In particular, Baylor et al. [3] conducted a neurophysiological experiment, which showed that the horizontal cells exert negative feedback on the photoreceptor cells [3].Then, electrical signals at the retinal ganglion cells are transmitted to the brain, passing through the optic nerve.However, the hierarchical mechanism of visual information processing in the human retina is not fully understood.To estimate the mechanism of visual information processing, electrophysio-logical experiments on animals and mathematical modeling and psychological experiments on humans have been actively conducted.
The mechanism of lateral inhibition is one of the most important functions in information processing in the early visual system.It is known that lateral inhibition is involved in understanding various early visual phenomena, such as visual adaptation and the Hermann grid illusion [42,47].In the field of mathematical modeling of the retina, a differential equation model assuming cell-cell interactions was proposed by [40] to understand the mechanism of lateral inhibition [16,40].In Peskin's model, negative feedback due to lateral inhibition induced by horizontal cells is expressed by the inhibition term in the following relation: where excitation is related to the light intensity.See "Appendix A" for a short review of Peskin model formulated as (4.1) and (4.2).In the first equation of (4.2) of Peskin's model, response, excitation, and inhibition in (1.1) correspond to u, u e and v, respectively.The properties of lateral inhibition are expressed by the second differential equation of (4.2).Subsequently, Krausz and Naka [22] and Tranchina et al. [49] proposed linear differential equation models of catfish retinal neurons [22,49], after which Peskin et al. [41] and Tranchin and Peskin (1988) proposed nonlinear differential equation models with nonlinear negative feedback [41,50].
Various electrophysiological experiments determined that retinal nerve cells have a radially two-dimensional symmetric receptive field on the retina.In particular, the retinal ganglion cells possess a Mexican-hat type receptive field, which consists of two parts: short-range activation and long-range inhibition [23].Based on the receptive field structure of the retinal ganglion cell, simple two-dimensional convolution models were proposed in [13,43].The following formula is a simple one-dimensional convolution model for the spatial variable x ∈ R: where I indicates an input signal, u indicates an output of the convolution model, and the integral kernel G * (x) indicates the Difference-of-Gaussian (DoG).The DoG is an example of a Mexican-hat-like function defined in Sect. 2. As shown in Fig. 1c, DoG is formed as where k ± are positive constants, and G ± (x) are Gaussian functions (Section 2.2 of [46]).In the convolution model (1.2), response, excitation, and inhibition in (1.1) correspond to u, R k + G + (x − y)I (y) dy, and R k − G − (x − y)I (y) dy, respectively.To describe the temporal response for an input signal I (t) depending on the time variable t ∈ R, a convolution model with temporal biphasic function G u (t) was proposed in [1,6].The definition of the temporal biphasic function is given in Sect. 2. Figure 1d shows an example of the temporal biphasic function G u (t) proposed in [1].As a convolution model with spatial variable x ∈ R and temporal variables t ∈ R, Adelson and Bergen (1985) proposed a nonseparable spatiotemporal oriented filter for motion perception in [1] as follows.They proposed a convolution model with (separable) integral kernel which comprises the Mexican-hat function G * (x) and temporal biphasic function G u (t) defined in the equation ( 1) of [1]; a number of (separable) integral kernels shifted in phase and time (Figs.9b-e in [1]) are summed to a spatiotemporal oriented filter (Fig. 9f in [1]).
After that, Lindeberg [26] proposed a mathematical framework constructed as the convolution model with a complicated spatiotemporal integral kernel described by the multiplication of a spatial affine Gaussian kernel and temporal smoothing kernel [26,27,29].Lindeberg's model can comprehensively explain the shapes of spatiotemporal receptive field structures of the retina, lateral geniculate nucleus (LGN), and primary visual cortex (V1) in the brain.
From the viewpoint of mathematical modeling by differential equations, for discussing the spatiotemporal integral kernel, Mahmoodi [31] proposed a nonlinear information processing model comprising several subunits constructed as (1.4) in the retina and explained spatiotemporal receptive field structures in LGN and V1 [31,32].One of the contributions of [31, p.139], was a condition given when the nonlinear information processing model behaved linearly.To derive subunits (1.4), partial differential equations were used.They derived a function G u (t) of the subunit (1.4) using v h (z, t) e N −1 e−1 (equation ( 9) in [32] and equation ( 5) in [33]), where N represents the number of spikes, which increases linearly with respect to the input potential v, i.e., N = kv, v h (z, t) is the fundamental solution for the following linear cable equation of the axon membrane's potential v: where c 1 is a positive constant and δ(t) is a delta function.If v is small, e N −1 e−1 is approximated by 1, and then the neural net- work behaves linearly.After that, an information processing model for the retina was designed on the basis of (1.4) as follows: (Fig. 15 in [31]), where g(x, s) is a series of Gaussian kernels, s = 0.3, ∂v h (z,t) ∂t is a temporal biphasic function, and is a Mexican-hat function (Fig. 1 in [28]).
The above mathematical models for understanding retinal information processing were applied to the development of image processing algorithms.The simple convolution model described by (1.2) has been applied to developments of image processing algorithms, such as edge detection [34].After that, to develop high-precision edge detection algorithms, partial differential equation models were proposed in [38,39].In [51], properties of some partial differential equations of image processing were explained comprehensively from a mathematical point of view, and some engineering applications of edge detection were shown.In Table 1 of [44], not only the correspondence between a diffusion equation and Gaussian, but also those between other partial differential or pseudodifferential equations of linear shift-invariant scalespaces and kernels of convolution models were summarised.In Section 1.3 of [45], a short overview about nonlinear anisotropic diffusion scale-spaces which can be used edge detections in images, were shown.Other partial differential equation models for image processing were also proposed in [9,18].
Recently, in the field of mathematical biology, novel mathematical methods for approximating time evolution equations with convolution integrals (1.2) to reaction-diffusion equations were proposed in [10,36,37].

Our Motivation
From the above short historical review of mathematical formulations for understanding the visual system, there are two major approaches: the differential equation model based on cellular interactions and the convolution model based on receptive field structures.Considering the above approaches, our study poses the following question: Which cellular mechanisms produce the spatiotemporal receptive field structure which consisting of the temporal biphasic function and Mexican-hat function?
Recently, to understand the derivation of Mexican-hattype receptive field structure, we proposed a hierarchical differential equation model that consists of the following stages based on the cellular interactions between the retinal nerve cells [46]: In this paper, we regard the whole structure of the retina as the one-dimensional space R, and differential equations of the Stages 1 and 2 are formulated as follows: where  1a shows a schematic diagram of our three-component model.
Next, we review our previous results.In [46], we considered the stationary problem for (1.5) and (1.6) and showed that the stationary solution u = u(x) can be represented as (1.2) and that the integral kernel G * (x) in (A.23) of [46]  The condition (1.7) implies that the lateral inhibition effect exists in our differential equations (1.5) and (1.6).In particular, we observed that the condition (1.7) is independent of the control parameter γ and investigated this property in Section 3.2 of [46].
In the subsequent study [20], based on our differential equations (1.5) and (1.6) for the light and dark inputs, we constructed a novel differential equation model for color inputs assuming the stage theory (trichromatic and opponentprocess theories), which is represented in the following form: The above form for the stage theory was proposed by Hall et al. [14], and the receptive field structures for the opponent colors were explained by Lindeberg [26].In [20], we studied the stationary problem of our model for color inputs in both cases of one-dimensional space R and two-dimensional space R 2 .Consequently, we showed that our model's output could be expressed in convolution integral form and derived the necessary and sufficient conditions for the appearance of the Mexican-hat function in our model's output.
According to the above review of our previous studies, we organized the derivations of the Mexican-hat function in our differential equations (1.5) and (1.6).However, the spatiotemporal integral kernel in our model has not been derived.In the neurophysiological experimental observations of the spatiotemporal receptive fields, a time lag was found before the Mexican-hat function appeared after excitation at x = 0. See Fig. 2 in [7], Fig. 1A in [2], and Fig. 4b in [30] for the bipolar cells and Figs.1C and D in [2] for the ganglion cells.Based on the above observations, we propose that at least the following two properties are satisfied (See Fig. 1b): (i) there exists t 1 > 0 such that K u (t 1 , x) > 0 holds for any x ∈ R, and (ii) there exists t 2 (> t 1 ) such that K u (t 2 , x) has a Mexicanhat shape.
Focusing on the above conditions (i) and (ii), in the convolution integral model (1.3) with (1.4), the integral kernel K u (t, x) has the Mexican-hat function for any t > 0. Thus, the condition (i) is not satisfied.Therefore, our motivation for this study is to derive the spatiotemporal integral kernel K u (t, x) of our model (1.5) and (1.6) and to investigate whether K u (t, x) satisfies the above two conditions (i) and (ii).

Our Contributions
The contributions of this study can be summarized as follows.
• In Sect.  .From these observations, we see that the time lag of the appearance of the Mexican-hat function in K u (t, x) depends on the time constant τ v .Therefore, K u (t, x) in our model satisfies conditions (i) and (ii) in Sect.1.1 under certain parameter conditions.• In Sect.3.1, we give numerical results for the appearance and non-appearance of a simple afterimage using the relation between the time constants shown in Corollary 1.
In Sect.3.2, focusing on an asymmetrical phenomenon called afterimage rotation, as investigated in [15], we introduce a nonlinear function τ v = φ(v) based on the asymmetrical temporal response of the inhibition exerted by the horizontal cells.We then give a theoretical prediction that the asymmetrical temporal response at the horizontal cell plays an important role in obtaining our visual impressions for the afterimages using Fig. 11.
The above results imply that a major factor in the construction of the temporal biphasic functions proposed in [1,6] is the temporal response speed of the horizontal cells.It is conjectured that the temporal response speed of the horizontal cells is a major factor for the construction of the spatiotemporal receptive field structure observed in [2,7,30].

Outline of this Paper
The remainder of this paper is organized as follows: In Sect.2, we first define temporal biphasic and monophasic functions and the Mexican-hat function.Sect.2.1 deals with the differential equations (1.5) and (1.6) with the input I = I (t) depending only on time t > 0 and gives a mathematical proof that the output u can be written in convolution integral form.In particular, the mathematical condition for classifying the integral kernel as the temporal biphasic or monophasic function is organized.Section 2.2 deals with (1.5) and (1.6) with input I = I (t, x) depending on time t and spatial variable x.The details of mathematical proofs for Main Results 1 and 2 are described in Appendix B.Moreover, Figs. 4 and 5 show that our spatiotemporal kernel possesses the temporal biphasic and Mexican-hat functions.Section 3 gives numerical results for obtaining visual impressions of the afterimages.Finally, Sect. 4 summarizes the paper.

Main Results
Here we give the definition of temporal biphasic and temporal monophasic functions according to [35,Fig. 3(a)].

Definition 1 A function G(t) (t > 0) is called a temporal biphasic function if the following two conditions are satisfied:
(i) There exists t = t 0 > 0 uniquely such that G(t 0 ) = 0 holds, (ii) G(t) > 0 holds for 0 < t < t 0 , and G(t) < 0 holds for t > t 0 .
We recall the definition of Mexican-hat functions according to [24].

Definition 3
We call a function f (r ) (r ∈ R) a Mexican-hat function if it satisfies the following conditions: (H 3) There exists a positive constant r such that f (r ) > 0 on an interval (−r , r ), and f The bottom (t = t 2 > t 1 ) of Fig. 1b shows a schematic diagram of a Mexican-hat function.
In order to make Theorem 2, Corollary 1, and Theorem 3 and these proofs easier to see, let us introduce We deal with (1.5) and (1.6), normalized at the parameters a, b and c as follows.Firstly, we divide the first equation of (1.6) by b, and replace τ u /b with τ u where note that we use the same notation τ u for the replacement τ u /b → τ u .Secondly, we apply the following replacements: Then (1.5) is rewritten as follows.
Moreover, we apply the following replacements: Then (1.6) is rewritten as follows.
In Sect.2.1, we consider that I = I (t) depends on only t > 0, and show that (2.1) and (2.2) can be written as and clarify when the integral kernel G u (t) to be the temporal biphasic or the temporal monophasic functions depending on two parameters C and R. In Sect.2.2, we consider an input signal I = I (t, x), depending on t and x.Our main result show that (1.5) and (1.6) can be written as an integral form with the kernel K u (t, x), and when t > 0 is sufficiently small, then the integral kernel K u (t, x) is positive for all x.The proof is written in "Appendix B".

Case I = I(t)
In this section, we consider that the input I = I (t) depends on the time t.Let us introduce the following notations: and can be written as (2.3), where where g 0 (t) and g 00 (t) are given by (2.18) and (2.19).
Proof of Proposition 1 Since I = I (t) depends on only t > 0, instead of (2.1) and (2.2), we deal with the following equations: (2.10) Firstly, we show that (2.9) can be written as the convolution integral  By integrating at the interval (0, t] the above equation, we obtain the convolution integral (2.11).
Secondly, we show that (2.10) can be written as the integral form where g 0 (t) and g 1 (t) are given by (2.18).Let us rewrite (2.10) as the following form: where the matrix A is given as follows.
When the matrix A has two distinct eigenvalues (2.5), we can solve Eq (2.14) by using the diagonalization for A as follows.Put then for λ + = λ − , there exists a inverse matrix such that By using (2.16), (2.14) is equivalent to (2.17) , then by using (2.6) and (2.15) we find that the solution of (2.17) is given by (2.18) Finally we show that (2.9) and (2.10) can be written as the convolution integral (2.3), where G 0 (t) and G u (t) are given by (2.8), where When substituting (2.11) for u e (s) in (2.13), the following integral appears.For the above integral, the following equation holds: (2.20) Note that (2.20) can be obtained as follows.By using (2.9), (2.12), the integration by parts, we have By substituting (2.11) for u e (s) in (2.13), using (2.20), we find that G 0 (t) and G u (t) of (2.3) are given by (2.8).It is easy to see that lim t→0 G u (t) = lim t→+∞ G u (t) = 0 holds.
Our results are the following Theorem 2 and Corollary 1.
Theorem 2 Assume that λ ± are real numbers, λ + = λ − , and η = ± .Then the followings hold.Remark 3 When λ ± are not real numbers, it is expected from numerical computations as shown in Fig. 3b that the shape of G u (t) is similar to damped oscillation.Therefore, G u (t) can be formed the temporal biphasic case, the temporal monophasic case, and the damped oscillation case, that is, this implies that our differential equation model possesses the integral kernel G u (t) with similar properties as the impulse response function which is well-known in the field of visual processing [6].

Remark 4
From Corollary 1 and Fig. 3a for profiles of the integral kernel G u (t), we can recognize that the time constant τ v plays an important role for the production of the biphasic function.

Proof of Theorem 2
Let us rewrite G u (t) as where . From (2.5), it is easily seen that λ − < λ + < 0 and the following holds. 1 From these and (2.7) we obtain Hence in the following we examine ( f 0 (t) − 1).By using the Taylor expansion and (2.7), we have From this and (2.21), we obtain Now we prove following inequality: Putting β 1 = + − η, and β 2 = − − η, then we find that For t > 0, we find that the following holds.

Proof of Corollary 1
It is obvious for (ii) and (iii), hence we show the proof only for (iii).It is easy to see that the inequality From this inequality and C > 1, we have R > C. Therefore, the condition of (iii) implies τ e < τ v and τ u < τ v .

Case I = I(t, x)
In this section, we consider that the input I = I (t, x) depends on the time t and the spatial variable x.We have the following Main Results 1 and 2 (accurate descriptions are Theorem 3 and Theorem 4 in Appendix B).

Main Result 2
When t > 0 is sufficiently small, then the integral kernel K u (t, x) of (2.26) is positive for all x.
The detail of the proofs of Main Results 1 and 2 are described in Appendix B.
When I (t, x) = δ(t)δ(x) (a strong bright stimulus is given for a moment in a controlled environment), the last term of (2.26) is equivalent to K u (t, x).Figs. 4 and 5 demonstrate spatiotemporal profiles of K u (t, x).In Figs.4a and 5a, K u (t, x) > 0 holds for small t, as shown in Main Result 2, and spatial profiles in K u (t, x) at t = 0.43 of Fig. 4b and at t = 0.3 of Fig. 5b are Mexican-hat type profiles corresponding to Fig. 1b.Also, temporal profiles in K u (t, x) at x = 0 of Figs.4c and 5c are formed as the shape of the temporal biphasic function.In particular, as shown in Fig. 5, in the case that the parameter τ v = Rτ u is small, the time lag for the appearance of the Mexican-hat function decreases, which Fig. 4 a A spatiotemporal profile of K u (t, x) with R = 10, where thick and dashed curves indicate contour lines at K u (t, x) = 0.005 and K u (t, x) = −0.005,respectively.A grayscale bar indicates the value of K u (t, x).b shows two spatial profiles in (a) at t = 0.4 (dashed line) and t = 0.43 (solid line), respectively.c is a temporal profile in (a) at x = 0.The other parameters are listed in Table 1 are compared with Fig. 4. Hence, it is suggested that the temporal response speed of the horizontal cell plays a vital role in obtaining the spatiotemporal receptive field observed in the field of neurophysiology.

Numerical Results
From the mathematical results shown in Sect.2, the temporal response speed of the component v plays an important role in the construction of the spatiotemporal integral kernel Fig. 5 a A spatiotemporal profile of K u (t, x) with R = 2, where thick and dashed curves indicate contour lines at K u (t, x) = 0.005 and K u (t, x) = −0.005,respectively.b shows two spatial profiles in (a) at t = 0.05 (dashed line) and t = 0.3 (solid line), respectively.c is a temporal profile in (a) at x = 0.The other parameters are the same as Fig. 4 with the temporal biphasic and the Mexican-hat functions.In this section, we demonstrate numerically that the temporal response speed of the component v is a major factor in the appearance of the afterimages.

Temporal Response of the Self-control Mechanism and Simple Afterimages
An adaptation is a transformed visual experience when we continue to observe a visual stimulus.Specifically, light visual stimuli are continuously observed for light adapta-tion, and that of the dark are continuously observed for dark adaptation.In general experiences, our eyes can adapt and be accustomed to higher than lower levels of environmental illumination.Light and dark adaptations may be a function of inhibitory feedback from horizontal cells to photoreceptors [25].When people adapt to a visual stimulus for a sufficiently long time, and then that visual stimulus disappears, the visual experience persists at the visual stimulus's location.This phenomenon is called an afterimage [19].It is thought to be explained by the spatiotemporal dynamics system of visual information processing within the retina.Barlow and Sparrock (1964) found that the longer the duration of dark adaptation, the lower the threshold of lightness perception and the lower the lightness of the perceived afterimage.This suggests that adaptation and afterimages may be caused in the same cells that process visual information.
Figure 6 shows a simple example of an afterimage.First, observers look closely at the center of a white disk on a black background shown in the left panel of Fig. 6 for a few seconds.Then, when the white circle is removed, and the black background is observed, the brightness within the region enclosed by the dotted circle in the right panel of Fig. 6 is perceived as darker than the outer region of the black background.Note that there is an individual difference for obtaining the above visual impression.
Next, we show numerical computations corresponding to the situation of an afterimage observed for the center of the white circular region in Fig. 6.Here, we consider the relation between the following input I (t) and output u(t) of our model (1.5) and (1.6) without the spatial variable x.
In our simulations, we use the implicit method with time discretization width dt = 0.001 and spatial discretization width dx = 0.01.Note that R is a free parameter.Figure 7a shows the output u(t) with R fixed as R = 10.0 > 1, that is, τ e = τ u = 0.1 < τ v = 1.0, which satisfies condition (ii) in Theorem 2. Hence, the integral kernel G u (t) is the temporal biphasic function, as shown in Fig. 3a.Because the time constant τ v is large, the effect of brightening the dark light stimulus has been operated after switching the input.Therefore, the value of u increases suddenly after switching from I 0 = 1.0 to I 1 = −1.0,and u decreases with time elapsed and has a stationary value.This is consistent with the afterimage observed in Fig. 6.However, when R is fixed as R = 0.1 < 1, as shown in Fig. 7b, G u (t) is the temporal monophasic function (see Fig. 3a), and the tendency of the afterimage cannot be confirmed in the profile of u(t).

Asymmetrical Temporal Response of the Self-control Mechanism and Motion Illusions by Afterimages
In recent years, motion illusions caused by afterimages have been studied mathematically and psychologically [15].Figure 8 shows an example stimulus of the Fraser-Wilcox illusion [11].The stimulus is perceived as clockwise and counterclockwise rotations in the peripheral vision.The perceived rotation depends on the pattern of brightness gradients [11,12], which is perceived in the direction of a change from black to white [17].When the Fraser-Wilcox stimulus is removed after a period of observation (adaptation), rotation is perceived in the opposite direction to that during adaptation.The rotation direction of this afterimage is known to be reversed when the background brightness is darkened.Hayashi et al. called such phenomena afterimage rotation [15].The previous study showed that the asymmetric phenomena in afterimage rotation could be explained by combining a convolution model with the temporal biphasic function and Weber's law, which represents psychological quantities.In [15], an output I modified by Weber's law is described by the function I = α log(I ) + 1, where I ∈ [e −1/α , 1] is an input, and α is positive.However, the visual mechanisms controlling the psychological quantities calculated by Weber's law are not easily comprehensible.Meanwhile, our model, which assumes spatiotemporal properties and interactions between the retinal nerve cells, has shown that the temporal response speed of the component representing the self-control mechanism is an important property for obtaining afterimages.Therefore, we focus on the temporal response speed of the horizontal cell.According to the neurophysiological experiment described in [5,8], the horizontal cells possess asymmetrical temporal response speed that the case of bright input is faster than the case of dark input.From this asymmetrical temporal response, we assume the following two properties for the time constant τ v : if the component v is positive (bright), then τ v is small, and if v is negative (dark), then τ v is large.Let us define the following nonlinear function τ v = φ(v) depending on v as a mathematical formulation with the above two properties: where σ 0 ≤ σ 1 and ρ are positive constants.Figure 9 shows a profile of the function φ(v) in (3.2) with parameters listed in Table 2.   1 and 2. In Fig. 10a, it is confirmed that u(t) converges to a stationary solution in a short time after switching from I (t) = 0.8 to I (t) = 1.0.However, u(t) converges to a stationary solution over a long time in Fig. 10b.Therefore, the introduction of the function τ v = φ(v) can represent asymmetrical temporal responses of the self-control mechanism.
Second, we try to explain the asymmetrical phenomenon of afterimage rotation using our model, assuming the asymmetrical temporal response of the component v. Here, we consider the simple input shown in Fig. 11a.It can be perceived as an afterimage in that if the brightness pattern disappears from the light background, the right part of a single gradation pattern is brighter than the right part and the afterimage moves in the left direction.In addition, if the pattern disappears from the dark background, the left part of a single gradation pattern is darker than the right part, and it can be perceived typically that the afterimage moves in the left direction.In particular, it has already known from [15]  that brightness for switching the left and right movements is slightly darker compared with the middle brightness.We consider a one-dimensional input I (t, x) in the bounded region [0, L] as follows: where N ∈ N, I * ∈ [−1, 1] is a constant indicating the brightness of the background, and I N , j is defined as follows: Fig. 11b shows an example of a profile I (t, x) = I N with N = 4, corresponding to the brightness pattern of Fig. 11a; and the parameters t 0 , t 1 , and L are listed in Table 3.
and x c = L/2 as the middle position of the interval N .We show numerical results in the interval 4 obtained as solutions of the (1.5) and (1.6) with the initial condition u e (0, x) = u(0, x) = v(0, x) = 0 for all x ∈ [0, L] and the periodic boundary condition.As a numerical scheme, we use the implicit method with time discretization width dt = 0.001 and spatial discretization width dx = 0.01.Let p (t) and p r (t) be peaks of profiles u(t, x) at the left and right sides in the interval N at time t ∈ [t 0 , t 1 ], where the left side of N is given as ) and the right is given as ).Let us define the minimum of the peaks p (t) on the left side [x 0 , x c ] as and the maximum of the peaks p r (t) on the right side [x c , x 1 ] as Let (t , x ) and (t r , x r ) be pairs of time and x-coordinate of P and P r , respectively.That is, P = p (t ) and P r = p r (t r ).We define an index for evaluating movements of afterimage profiles in our model as follows.
It is regarded that if A(I * ) is positive, the profile of the afterimage moves to the right, and if A(I * ) is negative, the profile moves to the left.Figure 12a shows a graphical image of the index A(I * ) with I * = 0, where the parameters of φ(v) are specified as σ 0 = 0.5, σ 1 = 1.0, and ρ = 10.0, and the input I is given as the case of N = 4.
Finally, in the input I defined by (3.3) with N = 4, we show numerical results regarding the dependency of the switched background brightness I * .Figure 12b shows profiles of the index A(I * ) for I * ∈ [−1, 1].Similar to the above result, when the parameters σ 0 and σ 1 have the same value, A(I * ) is a symmetrical profile.However, when σ 0 = 0.5 < σ 1 = 1.0,A(I * ) is asymmetrical and I * with A(I * ) = 0 is smaller than 0. This result is consistent with the previous study [15].Moreover, when σ 0 = 0.1 < σ 1 = 1.0,I * with A(I * ) = 0 is decreased more compared with the above case.Therefore, it can be considered that the parameters σ 0 and σ 1 represent individual differences in visual impressions.

Conclusion
We conclude that partial differential equations for visual information processing can be used to analytically show the shape of the spatiotemporal kernel, which represents a selfcontrol mechanism.
First, this paper provides the convolution integral forms of the output in our differential equation model with three components.In the case of the input I = I (t) that depends on only time t, Theorem 2 and Corollary 1 show that the construction of the temporal biphasic function derived in Proposition 1 depends on the relation of time constants of the three components.In the case of the input I = I (t, x), which depends on time t and spatial variable x, Main Result 2 in Sect.2.2 shows that the spatiotemporal integral kernel shown in Main Result 1 in Sect.2.2 does not have the Mexican-hat function for small t > 0.Moreover, it is shown in Figs. 4 and 5 that the shape of our spatiotemporal integral kernel is similar to that of the record data of the spatiotemporal receptive field at the retinal nerve cells [2,7,30].Therefore, we argue that one of the benefits of the differential equation model assuming cellular interactions is providing a conjecture for a main component responsible for plays the production of spatiotemporal receptive field structures.
Second, a time-dependent convolution integral model was used to identify the time constants important for the temporal response of retinal neurons.The convolution integral kernel with the temporal biphasic function was newly derived using a mathematical framework.
Third, this study successfully conducts the numerical simulation of afterimages, a visual phenomenon, using differential equations that take into account spatiotemporal response patterns.This finding indicates that afterimages can be predicted by differential equations assuming hierarchical visual information processing in the retina, without necessarily assuming the form of the kernel in the receptive field structure.We believe that this study is worthwhile because it lays the groundwork for investigating the quantitative nature of self-control mechanisms in retinal processing in conjunction with neurophysiological and psychological experiments.
where I indicates the input, u e indicates the excitation of a receptor, v indicates the inhibition of the receptor, and u is the response of the receptor.τ e and β are positive constants.
In our previous paper [20], it is shown that if I = I (x) depends only on x ∈ R, then the solution u of the stationary problem for (4.1) and (4.2) can be written as u(x) = R G * (x − y)I (y) dy with the following kernel: Furthermore, we showed that if I = I (t) depends only on t ∈ R, then the solution u of (4.1) and (4.2) for t ∈ R can be written as u(t) = G 0 (t) + t 0 G u (t − s)I (s)ds with the following kernel: where η = −  Taylor expansion of the integral kernel for t.When t is sufficiently small, then we can use the Taylor expansion of the integral kernel for t, however, when t is not small, then we can not use Taylor expansion of the integral kernel for t.Hence Theorem 4 holds true for only small t > 0.

B.1. Auxiliary lemmas
Now we show the definition of Fourier transformation and inverse Fourier transformation.For functions f 1 (t, x) and f 2 (t, ξ), let us denote the Fourier transformation of f 1 (t, x) with respect to x by f 1 (t, ξ), and the inverse Fourier transformation of f 2 (t, ξ) with respect to ξ by F −1 ( f 2 )(t, x), which are defined by It is well known that the following lemma holds, however, we recall the proof (see, [21]).
Lemma 1 Let τ > 0, d > 0 and consider the heat equation (4.5) with w(0, x) = w 0 (x).Then the solution of (4.5) is given by Proof Taking the Fourier transformation of (4.5), then we obtain By dividing this equation by dξ 2 , and putting τ = τ dξ 2 , I = I dξ 2 , then we have Hence, we have By taking the inverse Fourier of the w(t, ξ), we have Now we show the calculation of the second term of the righthand side (4.6), and we omit the calculation of the first term of the right-hand side (4.6).
Switching the order of integration, we obtain ), then by using complex function theory we have From Lemma 1, we find that the solution w of (4.9) is given by the following ) τu , (4.12) and K 1 is defined as follows:  Then we obtain (4.14) with the help of 1 2π R e ixξ dξ = δ(x).

B.2. Main Theorems
We show main theorems.The proof of theorems are shown in next subsection by using Lemmas 2 and 3.  where K u0 and K v0 are given by (4.15).
Theorem 4 There exists T > 0 such that K u (t, x) > 0 holds for x ∈ R and t ∈ (0, T ).

B.3. Proof of Theorem 3
We prove the former part of Theorem 3. Substituting (4.7) for u e (s, y) in (4.10), then, for example, one term becomes as follows.Switching the order of integration, we obtain

Stage 1 :
Reduction of the retinal image resolution, Stage 2: Lateral inhibition caused by non-local interaction, and Stage 3: Switching mechanism between lightness contrast and assimilation.
Fig. 1b.Next, we show numerical results with respect to the shape of K u (t, x) in Figs.4 and 5as two examples of τ v .We can see that the shape of K u (t, x) of Fig.4b t= 0.4 looks like the upper panel of Fig. 1b, and that of K u (t, x) of Fig. 4b t = 0.43 looks like the bottom panel of Fig. 1b.A similar situation holds for Fig. 5

is the temporal monophasic function. Corollary 1 Remark 2
If C and R satisfy the condition (ii), (iii), or (iv) in Theorem 2, then τ e < τ v and τ u < τ v hold.If R = 1, then λ ± are not real numbers.Therefore, we consider the case R = 1 in Theorem 2.

Fig. 3
Fig. 3 Profiles of the integral kernel G u (t) depending on the time constants τ e , τ u , τ v and the control parameter γ .The other parameters are fixed as 1.0.a A solid line indicates the case of temporal biphasic function with τ e = τ u = 0.1 < τ v = 1.0 and γ = 1.0.A dotted line indicates the case of temporal monophasic function with τ e = τ u = 0.1 > τ v = 0.01 and γ = 1.0.b The case of parameters specified as τ e = τ u = 0.1 < τ v = 1.0 and γ = 10.0

Fig. 6
Fig. 6 Simple example to demonstrate the effect of afterimage

Fig. 7 a
Fig. 7 a A profile of u(t), v(t), and I (t) with R = 10.0.b A profile of u(t), v(t), and I (t) with R = 0.1.In both figures, solid, dashed, and dotted lines indicate profiles of u(t), v(t), and I (t), respectively.Black arrows indicate the direction of the self-control effect by the component v

Fig. 8 aFig. 9 A
Fig. 8 a Brightness pattern that induces afterimage rotations.b The left figure has the same pattern as (a), and the right figure indicates the rotational directions of afterimage rotations via marked triangles

First, we give 2 Fig. 10 a
Fig. 10 a, b are profiles of u(t) and v(t) obtained in our model with (3.2), respectively.u(t), v(t), and I (t) are indicated by solid, dashed, and dotted lines, respectively.Thin arrows indicate durations to the steady state

Fig. 11 aTable 3
Fig. 11 a A two-dimensional brightness pattern corresponding to I 4 .b A profile of I 4 .Each arrow indicates the position of the interval j Table 3 Parameters in the input I (t, x) defined in (3.3) Parameter Value Description t 0 25.0 Switching time of I (t, x) t 1 50.0Finishing time of I (t, x) L 100.0 Length of pattern in I (t, x)
u and τ v indicate time constants, b indicates the speed that the output u becomes u e , γ indicates the control parameter, d v indicates the diffusion rate, β indicates the supply rate, and c indicates the decay rate.Figure τ 2.1, Proposition 1 shows that when the input I = I (t) depends only on time t, the solution u(t) for differential equations (1.5) and (1.6) is written as the convolution integral form with integral kernel G u (t).Theorem 2 clarifies when the integral kernel G u (t) becomes the temporal biphasic or monophasic function.Corollary 1 shows that if the time constant τ v is larger than both constants τ u and τ e , G u (t) is the temporal biphasic function.• In Sect.2.2, we show that the solution u(t, x) for model (1.5) and (1.6) is written in convolution integral form with the spatiotemporal integral kernel K u (t, x).Next, we show mathematically that T > 0 exists such that our spatiotemporal kernel satisfies K u (t, x) > 0 for x ∈ R and t ∈ (0, T ).This means that K u (t, x) in our model satisfies condition (i) in Sect.1.1 and the shape of K u (t, x) looks like the upper panel of