Abstract
Human brain response is the result of the overall ability of the brain in analyzing different internal and external stimuli and thus making the proper decisions. During the last decades scientists have discovered more about this phenomenon and proposed some models based on computational, biological, or neuropsychological methods. Despite some advances in studies related to this area of the brain research, there were fewer efforts which have been done on the mathematical modeling of the human brain response to external stimuli. This research is devoted to the modeling and prediction of the human EEG signal, as an alert state of overall human brain activity monitoring, upon receiving external stimuli, based on fractional diffusion equations. The results of this modeling show very good agreement with the real human EEG signal and thus this model can be used for many types of applications such as prediction of seizure onset in patient with epilepsy.
1. Introduction
Brain as the most complex organ in the human body controls all bodies’ actions/reactions by receiving different stimuli through the nervous system. Any stimulus stronger than the threshold stimulus is translated by the number of sensory neurons generating information about the stimulus and the frequency of the action potentials. After the action potential has been generated, it travels through the neural network to the brain. In various sections of the network and the brain, integration of the signals takes place. Different areas of the brain respond depending on the kind and location of stimuli. The brain sends out signals which generate the response mechanism.
During many years, numerous studies related to the brain response to external stimuli have been reported by scientists. Some researchers studied the brain response to different kinds of stimuli without proposing any model. In case of visual stimuli, we can mention the work done by Kaneoke et al. in analyzing the effect of the visual stimulus size on the human brain response using magnetoencephalography (MEG) [1]; see also [2, 3]. Other groups of researchers investigated the effect of auditory stimuli on the brain response. For instance, Will and Berg studied and compared the brain responses to periodic stimulations, silence, and random noise using electroencephalography (EEG) [4]; see also [5, 6]. Olfactory stimuli also were the main focus of some researchers. Sutani et al. investigated the brain response to pleasant and unpleasant olfactory stimuli using MEG signals. They found out that the MEG signals have recorded from frontal/prefrontal cortical areas of the brain has some differences in case of pleasant versus unpleasant stimuli [7]; see also [8, 9]. Different works have been reported on the investigation of the brain response to other kinds of stimulus such as emotional stimuli [10, 11] and pain stimuli [12, 13].
On the other hand, some scientists proposed some models of the human brain activity. On the microscopic level, the work done by Freeman in the modeling of the EEG arising from the olfactory bulb of animals during the perception of odors is noteworthy. He developed a set of nonlinear equations for this response which generates EEG like pattern [14, 15]. In another work Seetharaman et al. proposed a mathematical model for generation and propagation of action potential in a node of Ranvier and they called it as the phase-lagging model of single action potential [16]. When the microscopic models are extended to a macrolevel, then different methods are employed. Many of these models assume the cortical region to be a continuum. Liley et al. developed a set of nonlinear continuum field equations which described the macroscopic dynamics of neural activity in the cortical region [17]. These equations were used by Steyn-Ross et al. who introduced noise terms into them to give a set of stochastic partial differential equations (SPDEs). They also converted the equations governed by Liley et al. into linearized ODEs. This model could predict the substantial increase in low frequency power at the critical points of induction and emergence. They later used this model to study the electrical activity of an anaesthetized cortex [18–21]. Kramer et al. started with the equations given by Steyn-Ross and coworkers and neglected the spatial variation and the stochastic input. They believed that this gave rise to a set of ordinary differential equations (ODEs) for the modeling of the cortical activity. They showed that the results obtained from the SPDE model agree with clinical data in an approximate way [22], but they also stated that the spatial sampling of the cortex was poor because of inherent shortcomings in the equipment used. Kulish and Chan have suggested a novel method for the modeling of the brain response using fundamental laws of nature like energy conservation and the least action principle. The model equation obtained has been solved and the results show a good agreement with real EEGs [23].
Despite rapid advances in the studies related to the analysis of the human brain response, there has been less progress in the mathematical modeling of the human brain activity due to external stimuli. Yet, it seems that the contemporary level of developments in physics and mathematics makes establishing quantitative correlations between external physical stimuli and the brain responses to those stimuli.
This paper attempts to introduce a new mathematical model of the human brain response to external stimuli based on the fractional diffusion equation. At first we talk about the macroscopic level of brain organization and EEG signal as an alert state of human behavior monitoring. Then, by introducing fractional diffusion equations and considering the EEG signal as a fractal time series we model the EEG signal using mathematical equations. This model is then solved and after discussing different parameters in the model we provide some results and discuss about the model’s solution in details. Some concluding remarks are provided at the end.
2. Macroscopic Level of Brain Organization and the EEG Signal
In order to study the human neural activity, one can consider different levels of brain organization at different scales in time and space from a single neuron (microscopic level) to the whole brain organization (macroscopic level). In fact, the whole brain activity cannot be observed by measuring the activity of just a single neuron as far it is informative as it contributes to study the entire population of which it is a member.
In this research, we focus on the macroscopic level of brain organization. Macroscopic level of brain organization refers to the level of neural assemblies’ population in which each neural assembly interacts with other neural assemblies in close and distant cortical areas, exhibits spatial-temporal behaviour, and paints the human behaviour [24]. The functional behaviour of the brain is encoded in these spatial-temporal structures and can be extracted from the macroscopic quantities dynamics observed by EEG signal mostly [25].
During many years scientists have studied the human behavior by recording and analysis of the EEG signals from different areas of the brain. The EEG signal is the composition of different frequency bands (oscillatory activities-Alpha, Beta, …) which are structured coordinately (spatially, temporally). In fact, this signal has different characteristics that can be used in order to study the human brain response to external stimuli in our research. For instance, Figure 1 shows the EEG signal of a subject who received an external stimulus at s. In this figure the brain response to external stimulus can be seen as a sudden upward deflection at about s after the application of stimulus at s.
3. Fractional Diffusion Equation
Here, in order to develop our model we start with a simple equation:where
Equation (1) is the well-known diffusion equation where the coefficient is the diffusivity of the medium for the property . In fact, this equation arises in many descriptions of biological and physical phenomena, including Brownian motion [26], gradient driven chemical diffusion (with Fick’s law), and heat transfer (heat diffusion with Fourier’s Law). It is noteworthy that (1) is written in the Euclidean space.
On the other hand, the diffusion process can also been studied for fractals in fractal space. Fractals (such as random walk) may be defined as self-similar geometric objects whose scaling exponent (dimension) satisfies the Szpilrajn inequality:where is the scaling exponent (dimension) of the object and is its topological dimension, that is, Euclidean dimension of units from which the fractal object is built. In fact, Fractal and Euclidean geometries are conjugate approaches to the geometry of natural forms. Fractal geometry builds complex objects by applying simple processes to complex building blocks; Euclidean geometry uses simpler building blocks but frequently requires complex building processes.
Considering the diffusion process in the case of fractals, (1) is changed to (4) which is called the fractional diffusion equation [27]:where .
In (4), is the speed of propagation and is the Hurst exponent with the value within the range that brings predictability of signal into account. In fact, the Hurst exponent can be viewed as the probability of the diffusion process being persistent in a certain given direction. Note that the case , which corresponds to a nonfractal diffusion process, leads to the well-known classical equation. Observe also that if , (4) degenerates into the Poisson equation, ; that is, there is no preferred direction of random walk in this case, while the case leads to the wave equation, .
It is now necessary to make a very important remark. It is possible to consider a Brownian motion type process as a process which takes place in an Euclidean space (see the right hand side of (4)), considering the temporal dimension, , of the process as fractal time which, for the same diffusion coefficient, , either slows down (in case of ) or speeds up (in case of ) the process in question. This can be described with the generalized diffusion equation (4) in which the time coordinate appears as a fractal quantity. The generalized diffusion equation is a fractional PDE of order with respect to time.
Taking all these conjectures into account, (4) is a generalized form for describing the diffusion process when the time becomes fractal.
4. Fractional Diffusion Model of EEG Signal
In this section, it is aimed to show that EEG signal as a fractal time series (but not stochastic), which represents a transient record of a random walk process, can be modeled by the solution of the fractional diffusion equation:where the term , which stands for the brain response, is the voltage fluctuations resulting from ionic current flows within the neurons of the brain. The impulse propagation speed within the neural network, which is a finite quantity, is represented by . The term , the effective diffusion coefficient, as the property of neural tissue, is related to neuron’s resistance to the electrical impulse as it travels over the nerve. The term , as a time variable parameter, corresponds to the Hurst exponent that brings the predictability of signal into account.
The direction of deflection at each moment in the signal can be studied by computing the Hurst exponent (Figure 2). The Hurst exponent is an indicator of the long term memory of the process generating the signal and thus it is the measure of the predictability of the signal.
As it was mentioned before the Hurst exponent can have any value between 0 and 1, where the value that it gains in each moment determines the behavior of the next deflection in the signal.
Firstly, if the Hurst exponent has a value between 0 and 0.5, it means that the process is antipersistent; that is, the trend of the value of the process at the next instant will be opposite to the trend in the previous instant. Secondly, a value of between 0.5 and 1 means that the process is persistent; that is, the trend of the value of the process at the next instant will be the same as the trend in the previous instant. Finally, if , the process is considered to be truly random (e.g., Brownian motion). It means that there is absolutely no correlation between any values of the process.
One of the interesting points about the fractional diffusion model (5) is that it accounts for a finite time lag (reaction time) between any given disturbance (stimulus) and the brain response to it (human action/reaction) based on the assumption that no instantaneous propagation of information is possible within the brain. This effect is considered during the derivation of this model, but it is substituted by , according to the equation:Equation (6) is the formula for relaxation time in hyperbolic reaction diffusion equations [28].
In this research we consider the brain in an informational space where is the spatial dimension of this space. In fact, each two neurons in informational space, which have informational interlink, make an informational channel between themselves and the length of this channel is called the informational distance where
It is noteworthy that informational distance is not same as the concept of spatial distance. Two neurons, which are close together, may not exchange information which means . On the other hand, two neurons which are far apart may be closely interlinked to each other and exchange a lot of information which makes very small. Thus, has value in the range [29].
The exchange of information between two neurons in the brain happens when there is a potential difference along the channel. In fact, the potential difference between two neurons causes an information flux and the exchange of information from the neuron with higher informational potential value to the one with the lower value continues till the gradient becomes zero.
The fractional diffusion equation is valid for the time scale:Also, regarding (5),is the initial condition at . Also, the solution of the model is required to be bounded as . Otherwise, the conservation of energy principle would be violated. In other words,is the boundary condition in the case of infinite domain.
In order to solve the above fractional diffusion model the method proposed by Oldham and Spanier is employed here [30].
Upon introducing the excess value , so the initial condition imposed on is .
In order to solve (5) we apply the Laplace transform with respect to time, ; then we havewhere is the Laplace transform variable and denotes the Laplace transform of the excess value .
Equation (11) is a second-order ordinary differential equation (ODE), where the general solution is calculated aswhere and are two arbitrary functions. However, since the solution is to be bounded for all , the second arbitrary function, , must be identically zero, so the solution (12) is changed toUpon differentiating (13) with respect to ,After comparing (13) and (14), can be eliminated; then it can be written asBy taking the inverse Laplace transform of (15) and restoring the original variables, then we havewhich is written in terms of fractional derivative of order with respect to .
Using the definition of fractional derivative [30], namely,where is the Gamma function, and noticing the Fick’s law,where represents the flux, and (16) can be written asSince the function represents the flux it can be equated with the external influence acting on the system (external stimulus).
Equation (19) provides the relationship between the nonequilibrium value, , and the external influence acting on the system, . This equation is valid for every location within the domain (including the boundary) at every moment.
Since is a nonnegative parameter, it follows from (19) that the value of , on the average, increases with the time according to the power law as , provided of course that the fluctuations are small in comparison with the averaged influence. Note that, for , and (19) yields the well-known diffusion (random walk) growth given by .
An external influence can be modeled by a Gaussian pulse; that is,where denotes the moment, at which the Gaussian pulse reaches its maximal value , whereas is the standard deviation of the Gaussian pulse.
In case of different types of stimuli, depending on the size and duration, the parameters in (20) will have different values. In the case of many concurrently external stimuli, a series of Gaussian pulses is considered.
It is noteworthy that (20) is not a compulsory form of external stimulus where other formulas might be used, provided that they meet the principle demands discussed in this research. This mathematical form is chosen because its capability was examined in the modeling of generation of action potential in a neuron [16].
In the next section we provide a formulation for computing the effective diffusion coefficient at each moment that can be used in (19) in order to compute the brain response to external stimuli.
5. A Phase-Lagging Diffusion Based Model of the Diffusion Coefficient
In order to analyze the behavior of fractal time series, the value and the direction of each fluctuation in the signal should be analyzed. The direction of each fluctuation can be found by computing the value of the Hurst exponent at the previous point of the signal.
In order to know about the value of the signal, we should know about different parameters which are appeared in the fractional diffusion model. By looking at (19) and replacing the value of , we have
The value of the Hurst exponent which is used in (21) can be computed using MATLAB based on Rescaled Range Analysis method [31, 32]. Then, in order to compute the value of the signal in (21) two parameters and should be known. In order to do this, first, a relationship between and can be made, and then by using the formula for relaxation time, (6), can be replaced by and .
In order to make a relation between and , the phase-lagging model of action potential can be used [16]:
In (22) the diffusivity term is related to the resistance of the neuron to the electrical impulse. is the property of the neural tissue and will dampen the impulse as it travels over the nerve. The term is the reaction time and is the zero-order modified Bessel function.
Considering the phase-lagging model of action potential in the whole brain scale then, this model and the fractional diffusion model explain the same phenomenon. Thus, the mathematical equations that belong to these two models can be equal.
By writing this equality,where is the diffusion coefficient in the phase lagging model of action potential and is the diffusion coefficient in the fractional diffusion model.
Considering and removing from both sides of the equation,Then, we can writeBecause (25) is valid for all values of then we can remove the integral:Thus,Diving both sides of (27) by and introducing a new variable we have
Considering a constant value for , in order to study the variation of versus , Figures 3 and 4 are provided. It is noteworthy in order to analyze the value of the signal that it is only needed to concentrate on the calculations of and, accordingly, in or because as it is known, having the probability of on one side yields the probability of in the reverse side, but with the same value of fluctuation. In this research the span is taken as the reference. So, for instance, in order to compute the value of the signal in a point with it is needed just to compute the value of the signal with .
As it can be seen in Figures 3 and 4, first the value of dimensionless diffusivity increases but after that as time increases (), dimensionless diffusivity tends a constant value. In order to describe this behavior, (28) can be analyzed when : is computed asBy substituting (30) into (29),Considering (31) in case of which is the same as the trends observed in Figure 3.
By considering (31) in case of we have Thus,which is the same as the trend observed in Figure 4 for values of where . So, in both cases ( and ) the dimensionless diffusivity tends to a constant value as time increase.
So, in cases of Figures 3 and 4, before the maximum point of graph, as time goes on, the value of the dimensionless diffusivity increases, but after passing the maximum point, the dimensionless diffusivity shows the opposite behavior and, after some time, the dimensionless diffusivity tends to a constant value which is 0 and 1 in the cases of and , respectively.
As has been considered as a time dependent parameter and as a constant, thus from previous discussion it can be concluded that when human senses a stimulus, the diffusion of this stimulus to human brain increases but after some time, the diffusion of the stimulus decreases. In fact, these results validate and show the importance of (28) for computing the value of .
Thus, by computing the value of at each time moment and substituting its value in the fractional diffusion model the value of the signal at each time moment can be computed.
6. Result and Discussion
In this section using the fractional diffusion model the human brain response to a visual external stimulus is modeled and compared with the real EEG signal.
6.1. Subjects
In this research the experiments were carried out on 6 voluntary healthy students (21–24 years old, 3 males and 3 females). Prior to the experiment, each subject was examined and interviewed by a physician to ensure that no neurological deficit, pain condition, or medication affects the EEG. All subjects were right handed. Informed consent was obtained from each subject after the nature of the study was fully explained.
All procedures were approved by the Internal Review Board of the University and the approval for experimentation involving human subjects was issued. It is noteworthy that the identity of all subjects remains confidential.
6.2. EEG Recording
The EEG data used in this research were collected using Mindset 24 device, a 24-channel topographic neuromapping instrument, which can measure 24 channels of data with the sampling frequency of 256 Hz.
In an electrically shielded, acoustically isolated, and dimly illuminated room a visual stimulus applied on subjects. It should be mentioned that it is endeavoured to insulate the subjects from all other external stimuli. This ensures that the response measured in the EEG signals is primarily due to the stimulus applied. In the experiment subjects were watching a checkerboard pattern (see Figure 5) on the monitor of a computer from the distance of 130 cm.
The stimulus was the checker reversal. After one second the reversed pattern (see Figure 6) was displayed and then the original pattern was displayed again. Thus, the stimulus was applied at s. The interstimulus interval of 5 s was chosen in these experiments.
Mindmeld 24 software is used for the collection of data using Mindset 24 machine. The software gives data in the form of .bin files which can be processed to give text files (.txt) that are required for further processing.
Although EEG data are recorded from 24 electrodes, in this research the analysis is done on the data governed from the left occipital (O1) electrode (near to the location of the visual primary sensory area). This electrode was chosen based on the nearest place to the visual sensory area which shows the strongest response that can be seen in the signal recorded from this electrode compared to other electrodes. The electrode impedance was kept lower than 5 KΩ.
A bipolar electrooculogram (EOG, vertical and horizontal) was recorded for off-line artifact rejection. After bandpass filtering in the range of 0.1–70 Hz, 2 seconds of data (256 data before stimulation and 256 data after stimulation) was saved. It means that there are 256 values of voltage collected every second. Prestimulation is defined as the status before the application of the stimulus. On the other hand, poststimulation is defined as the status after the application of the stimulus. As it was mentioned previously the stimulus was applied at s.
It is noteworthy that choosing higher sampling rate will result in more pre- and poststimulation data. Having more data will results in more precise computation of Hurst Exponent, diffusion coefficient, and accordingly the signal and its parameters such as response initiation time.
In the first week 40 trials were collected from each subject in one day. The data collections were repeated after a week for each subject in order to examine the reproducibility of the results from experiments. By repeating the experiments in the second week totally 80 trials were collected. After visual inspection of data collected from each subject and rejection of trials with artifacts, 40 trials free of artifacts were selected for future analysis. It is noteworthy mentioning that physician monitored the subjects during all experiments.
6.3. Data Analysis
A set of codes was written in MATLAB software in order to compute all required parameters which were discussed before.
As the recorded data were noisy, the EEG signals were filtered using the Wavelet toolbox in MATLAB and then were processed by the methodology discussed in this section.
The value of can be read from the record of the signal at the moment the stimulus is applied to the subject, s. The initial value of is computed for 1 second of the recorded data before the application of the stimulus to the subject. In order to compute the Hurst exponent, as it was mentioned previously the Rescaled Range Analysis method is employed, which is widely used by statisticians.
It is required to compute the Hurst exponent value in each moment in order to analyze the generated signal. At the first step, the program computes the Hurst exponent for the recorded EEG signal and the predicted signal and generates two time series in one figure.
In each moment, the program computes using (28). Also, as it was mentioned previously that, for all analysis performed here, a single stimulus is considered and thus a single Gaussian pulse is modeled using (20). By substituting the required parameters in (35), the program computes the value of the signal in each moment and then plots the modeled signal in a figure together with the recorded EEG signal (after stimulation):It is noteworthy that (35) is governed by substituting into (19).
The values of some required parameters are listed in Table 1.
As it was mentioned previously, one second of EEG data was recorded, and using these data as the reference, one second of EEG time series is predicted. After that, the modeled signal is analyzed in terms of the initiation time for the response, the response duration, and the peak to peak voltage. In this research in order to see the response fluctuations clearly, the real or the modelled signals after the poststimulation are averaged in each case. It means that for each subject 40 selected trials are used as the input to the model and, accordingly, the grand average of all modeled signals (40 signals for each subject) is presented. Then, the generated plots for the real or the predicted signals are the result of averaging over all selected trails in each case.
Here it should be mentioned that in the analyses of the EEG signals plots the initiation time for the response and the duration are chosen based on the literature notes which consider the major positive or negative pole. For instance, the fluctuations which have voltage in the range of 5 to 10 V or −5 to −10 V are related with the response to the stimulus. Thus, for instance, in order to have a feature of the response duration, the time span between the first peak and the last peak which have the voltage values within one of these two ranges is considered.
The grand average of the recorded EEG signals and the grand average of the predicted signals using the fractional diffusion model over all selected trials for 1 second post-stimulation are shown in Figure 7.
(a) Subject 1
(b) Subject 2
(c) Subject 3
(d) Subject 4
(e) Subject 5
(f) Subject 6
As it can be seen for different subjects in Figure 7, the recorded signal (black solid line) and the predicted signal (red dashed line) show the similar behavior. In both cases the brain response to the stimulus starts with a positive peak () after the application of the stimulus to the subject where its amplitude goes further than 5 V. This response causes the signal’s voltage to fluctuate in a bigger span. Following the positive peak a negative rebound () can be seen in the plot. In fact, the response to the stimulus terminates at this negative peak, after which the brain goes back to its normal status during rest, without any big deflection in the signal.
As an example, for subject 1 in cases of the real and predicted signals the response to the stimulus starts with a positive peak () at about 118 ms and about 127 ms, respectively, after the application of the stimulus to the subject. This response causes the signal’s voltage to fluctuate in a bigger span. Following the positive peak a negative rebound () at about s and s can be seen in the plot in cases of the real and predicted signals, respectively. In fact, the response to the stimulus terminates at this point, after which the brain goes back to its normal status during rest, without any big deflection in the signal.
The values of peak to peak voltage, the initiation time for the response (the initial peak of the response), and the response duration (the time difference between the first and the last peaks within the response duration) for the recorded EEG signal and the predicted signal in case of different subjects are provided in Table 2.
As it can be seen in Table 2 for all subjects the predicted initiation time for the response, the response duration, and the peak to peak voltage have very close values with their related values in the recorded EEG signals. Thus, it can be said that the predicted signal resembles the real EEG signal within the response duration in the cases of the initiation time for the response, the peak to peak voltage, and the response duration. Moreover, in order to study the uncertainty and predictability of the model’s solution, the Hurst exponent variations for the recorded EEG signals and the predicted signals over all selected trials for 1 second after stimulation are shown in Figure 8.
(a) Subject 1
(b) Subject 2
(c) Subject 3
(d) Subject 4
(e) Subject 5
(f) Subject 6
The high correlation between the values of the real signals and also the predicted signals can be realized by looking at the values of the Hurst exponents. For instance, in the case of subject 1, the value of the Hurst exponent is distributed between 0.900 and 0.943 for the recorded EEG signal (black solid line) and between 0.900 and 0.952 for the predicted signals (red dashed line). Thus, the low uncertainty of the prediction can be confirmed, and it can be said the signal is predicted well, because the Hurst exponent values are not close to , which stand for a truly random process. This behavior can be seen in the Hurst exponent plots for all subjects.
Also, as it can be seen in the Hurst exponent plot, the value of the Hurst exponent in the case of the real EEG signal and the predicted signal experiences a sudden upward deflection. For instance, in the case of subject 1, the value of the Hurst exponent is decreasing in the time span of s to about s and s to about s, respectively, for real EEG signal and the predicted signal; after that, a sudden upward deflection can be seen, which stands for experiencing the visual stimulus, and again the trend shows the same behavior. The overall decreasing behavior stands for the phenomenon that when a longer time span is considered, the less the human brain “remembers” its initial state. The same behavior can be seen in other Hurst exponent plots.
By analyzing the behaviors which have been seen in Figures 7 and 8 it can be said that on one hand the uncertainty of the prediction was low and on the other hand, the accuracy of the prediction was very good as the predicted signal resembles the real EEG signal.
All the analyses which have been done in this research show that EEG signals can be modeled by the solution of fractional partial differential equations and, thus, the behavior of system modeled by means of such equations can, in principle, not only be predicted but also quantifies.
7. Conclusion
In this paper we introduced a new mathematical model which quantifies the human brain response to external stimuli. We developed this model by applying the fractional diffusion equation to human EEG signals. The model generates a multifractal time series which shows a quantitative concurrence with the real EEG signals. Using this model we successfully predicted the EEG signals of different subjects upon receiving a visual stimulus. This model shall be further applied in case of different external stimuli where the results can be verified against the real EEG signal which means the prediction of the human behavior by forecasting the EEG signal. On the other hand this model also can be employed in order to predict different abnormal brain activities, such as epileptic seizure, by at least some seconds before the time of occurrence. If so, a seizure warning and the expected time of this epilepsy occurrence can be generated, leading to the future monitoring of this disease.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.