A computational model unifies apparently contradictory findings concerning phantom pain.

Amputation often leads to painful phantom sensations, whose pathogenesis is still unclear. Supported by experimental findings, an explanatory model has been proposed that identifies maladaptive reorganization of the primary somatosensory cortex (S1) as a cause of phantom pain. However, it was recently found that BOLD activity during voluntary movements of the phantom positively correlates with phantom pain rating, giving rise to a model of persistent representation. In the present study, we develop a physiologically realistic, computational model to resolve the conflicting findings. Simulations yielded that both the amount of reorganization and the level of cortical activity during phantom movements were enhanced in a scenario with strong phantom pain as compared to a scenario with weak phantom pain. These results suggest that phantom pain, maladaptive reorganization, and persistent representation may all be caused by the same underlying mechanism, which is driven by an abnormally enhanced spontaneous activity of deafferented nociceptive channels.

Phantom experiences are vivid sensations of a body part that was lost after an accident or surgery. These experiences are very common among amputees, and 80% of them report intensely painful sensations1 which are commonly comprised under the term phantom limb pain (PLP). To explain phantom limb pain, a model of maladaptive reorganization has found considerable attention234. The model is motivated by the finding that the amount of reorganization of the somatotopic map in the somatosensory cortex S1, the somatosensory homunculus, strongly correlates with the phantom pain rating of patients. Using functional magnetic resonance imaging (fMRI) or magneto-encephalography (MEG), the reorganization is usually assessed as the spatial invasion of adjacent cortical representations in the somatotopic map into the representation of the missing limb. The model of maladaptive reorganization has sparked the development of novel therapies of phantom pain5678910.

Recently this model has been questioned, since fMRI measurements during executed movements of the phantom limb – which are to be distinguished from merely imagined movements11 – revealed that the cortical representation of the missing limb was still preserved and showed an even more pronounced activation than in patients without phantom experiences12. Moreover, the measured activity of the cortical representation of the phantom limb was positively correlated with the phantom pain rating of the patient. These new findings seem to contradict the model of maladaptive reorganization and motivated the authors of the study to propose an alternative model of persistent representation which is currently debated1213.

The aim of the present study is to investigate the conflicting findings concerning the cortical representations in phantom pain by means of two variations of organizing the cortical map using a computational model.

The model (Figure 1) is based on the following physiologically plausible assumptions:

Figure 1

Schematic representation of two variations of the computational model: (a) Model variation A involving an integrated cortical map receiving input from both tactile and nociceptive channels.

(b) Model variation B involving a split cortical map receiving separate input from tactile and nociceptive channels. The tactile and nociceptive channels (red and blue lines) sample the entire skin surface. For simplicity, a homogeneous distribution of receptors was assumed, and only the fingers of one hand were modeled. Each channel has a receptive field on the skin with a specific threshold and gain for either tactile or nociceptive stimuli. There are three types of gate per channel: peripheral gates, spinal gates, and central gates. The peripheral gates control the receptor sensitivity, the spinal gates control the afferent flow through the dorsal horn, and the central gates control the input into a self-organizing cortical map. Each channel is independently afflicted with discrete neuronal noise (DNN) and spontaneous coherent activity (SCA). For the simulations, the only modified parameters after amputation correspond to the channels originating in the amputated finger. These are 1) the spinal gate threshold, 2) the SCA strength, and 3) the central gate threshold. See Table 1 for an overview over the model parameters on the different conditions.

A1) The somatosensory cortex is a neural network that dynamically organizes itself according to incoming neural activity. In the computational model, the self-organizing network is implemented as a Kohonen map14.

A2) There is spontaneous activity in the sensory system, which is under normal circumstances weaker than the regular activity induced by actual stimulation, and which is abnormally increased in those parts affected by deafferentation. In the computational model, spontaneous activity in the sensory pathway is implemented in a twofold manner: There is discrete neuronal noise (DNN) in the form of randomly occurring discrete events with randomly varying amplitude, and there is spontaneous coherent activity (SCA) in the form of randomly occurring activation events that resemble those events caused by actual stimulation.

A3) Movement execution of, and attention to, a (phantom) limb activates networks in the brain that have stored sensory experiences of the (phantom) limb. In the computational model, voluntary movement of the (phantom) limb would enhance the SCA.

A4) The somatosensory afferent input to the cortex is regulated by a neural mechanism analogous to the gate control theory of Melzack and Wall151617. Specifically, the hypothesized mechanism compensates for the long-term input strength, so that a long-term increase/decrease of input strength would eventually lead to an increase/decrease of the gating threshold. In the computational model, the regulation is implemented as a linear, saturated gate, the central gate. There are two more, functionally analogous, gates implemented in the model, the peripheral gate and the spinal gate, which control the receptor sensitivity and the afferent flow through the dorsal horn of the spinal cord, respectively. Also the spinal gate responds to a long-term increase/decrease of input with an increase/decrease of its threshold.

For both model variations, the map would organize itself according to topological information only (where-information), and would not take into account the modality (touch or nociception) of the signal (what-information). Therefore, the map encodes only where-information and no what-information. For model variation A involving the integrated cortical map receiving input from both tactile and nociceptive channels, only the summed activity of tactile and nociceptive channels influences the self-organization process. For the model variation B involving the split cortical map receiving separate input from tactile and nociceptive channels, each sub-map is organized according to either tactile or nociceptive input, thus reflecting the proportional activity from the topological regions of one specific modality.

Results

The simulations of the computational model have been run on different sets of parameter values for each of the different model variations (integrated and split map). These sets are denoted as PRE, NOPAIN, and PAIN, and they are taken to represent different physiological conditions (Table 1). The internal state of the model resulting from a simulation on the PRE condition served as the starting point for the simulations on the NOPAIN and PAIN conditions, respectively.

Table 1

Values of the model parameters on the three conditions PRE, NOPAIN, and PAIN, for the tactile and nociceptive channels (“Tact” and “Noci”). On the PRE condition, all channels share the same values (“All fingers”). On the NOPAIN and PAIN conditions, only those channels originating in the amputated finger are modified with respect to the PRE condition, and only those are indicated here (“Amputated finger”). Blank fields correspond to values that are unchanged with respect to the PRE condition. Values in brackets correspond to the values taken in the probing phase. Abbreviations are as follows: Tact = tactile channels; Noci = nociceptive channels; stim = stimulation; amp = amplitude; dur = duration; DNN = discrete neuronal noise; SCA = spontaneous coherent activity; p-Gate = peripheral gate; s-Gate = spinal gate; c-Gate = central gate

	PRE		NOPAIN		PAIN
	All fingers		Amputated finger		Amputated finger
Parameter	Tact	Noci	Tact	Noci	Tact	Noci
stim rate	0.2	0.01	0	0	0	0
DNN rate	2	2
SCA rate	0.2 [1.0]	0.01 [0.05]				0.05 [0.25]
stim amp	1	1
DNN amp	0.05	0.05
SCA amp	0.05 [0.25]	0.05 [0.25]				0.25 [1.0]
stim dur	0.1	0.1
SCA dur	0.1	0.1
p-Gate thresh	0.1	0.1
s-Gate thresh	0.1	0.1	0.025	0.025	0.025	0.025
c-Gate thresh	0.1	0.1	0.025	0.025	0.15	0.15
p-Gate gain	1.234	1.234
s-Gate gain	1.234	1.234
c-Gate gain	1.234	1.234

The PRE condition is taken to represent the situation of a healthy subject before amputation. Spontaneous activity in the somatosensory channels is at baseline level and the gate thresholds are reasonably low, so that only very few spontaneous activation events (DNN and SCA) are strong enough to pass the spinal and central gates and reach higher regions of the brain including the somatosensory cortex. In line with assumption A1, the somatosensory map organizes itself into regions that correspond to parts of the skin in a somatotopic manner (Figure 2, top row). According to assumption A3, the execution of a movement of the middle finger enhances the SCA in the channels corresponding to the moved finger. Consequently, some SCA events are just strong enough to pass the central gate and contribute to a weak central activation of the cortical region corresponding to the moved finger (Figure 3, top row).

Figure 2

Representative example of the cortical map in model variation A, after the training phase with the parameters from the three conditions PRE (before amputation), NOPAIN (after amputation of the middle finger, without SCA enhancement), and PAIN (after amputation, followed by strong SCA enhancement).

The example is taken from thirty simulations of model variation A involving the integrated cortical map. Colored units in the map (left panel) are associated with receptors in equally colorized regions on the hand (right panel). Black units in the map are not associated with any of the receptors. Simulations of model variation B lead to similar maps with the exception that the proportions of the finger representations in the nociceptive map evolve differently from those in the tactile map, as numerically shown in Figure 5.

Figure 3

The same representative simulation of model variation A as in Figure 2, showing the activity of the integrative cortical map, summed over the probing phase (movement of the existing or phantom middle finger) after training with the parameters according to the conditions PRE, NOPAIN, and PAIN.

The representations of the activations on the hand are shown as colored dots on the hand. The color encodes the strength of the activity in arbitrary units (accumulated normalized firing rate). Simulations of model variation B lead to activities represented in Figure 4.

The NOPAIN condition represents the situation of the same subject as in the PRE condition after amputation of the middle finger. According to assumption A4, the loss of afferent input to the affected sensory channels causes a decrease of the threshold of the spinal and central gates. Due to the decrease of the thresholds, more spontaneous activity events pass the spinal and central gates and enter the somatosensory cortex. Consequently, the representation of the missing finger is preserved in spite of the absence of external sensory input (Figure 2, middle row). As the spontaneous activity events passing the central gate remain weak, the amount of central activation caused by the spontaneous activity is small, even during the execution of phantom movements (Figure 3, middle row).

The PAIN condition represents the situation of the same subject as in the PRE condition after amputation of the middle finger and with permanently increased spontaneous activity in the affected nociceptive channels. According to assumption A4, the increased spontaneous activity entering the central gate would lead to an increase of the gating threshold, so that only few but very strong SCA events would pass the gate and activate the corresponding region of the somatosensory cortex. The net result is a cortical map with a smaller, but preserved, cortical representation of the phantom (Figure 2, bottom row). As the spontaneous activity events passing the central gate are very strong, the cortical representation of the missing finger is strongly activated during the execution of phantom movements (Figure 3, bottom row).

A statistical comparison of the simulation results for model variations A and B, and on the individual conditions PRE, NOPAIN, and PAIN, yields:

In the absence of any stimulation, the central activity from all sensory channels was zero on the PRE condition (Figure 4a). On the NOPAIN and PAIN conditions, the central activity from the nociceptive channels of the amputated finger were nonzero, and they were significantly higher on the PAIN condition than on the NOPAIN condition (Figure 4a, dark gray bars). On all conditions, the central activity remained zero for sensory channels of the non-amputated fingers (data not shown).

Figure 4

Central tactile (white), nociceptive (dark gray), and total (light gray) activity from sensory channels originating in the (amputated) middle finger.

(a) Central activity accumulated over the resting phase (no stimulation), and measured in terms of the normalized firing rate (arbitrary units). The values depicted are median values, so the total central activity (light gray) is the median of the sum of both tactile and nociceptive activities. The activities for the remaining fingers (not depicted) where zero on all conditions. (b) Central activity accumulated over the probing phase (during voluntary phantom movement), also measured in terms of the normalized firing rate (arbitrary units). Error bars indicate the (25,75)% quantiles, braces indicate significant differences (***: p < 0.001, Wilcoxon rank sum test, statistically corrected). The shown values result from simulations of model variation A. Model variation B is mathematically equivalent with respect to the depicted values.

For model variation A, there was significantly more cortical reorganization on the PAIN condition than on the NOPAIN condition (Figure 5a). For model variation B, there was significantly more reorganization in the tactile part of the cortical map, on the PAIN condition than on the NOPAIN condition (Figure 5b, white bars). On the other hand, there was significantly less reorganization in the nociceptive part of the cortical map, on the PAIN condition than on the NOPAIN condition (Figure 5b, dark gray bars). Also, the nociceptive representation of the amputated limb significantly shrinked on the NOPAIN condition as compared to the PRE condition (Figure 5b, dark gray bar on condition NOPAIN is negative), while it did not significantly change on the PAIN condition (Figure 5b, dark gray bar on condition PAIN is not significantly different from zero).

Figure 5

Cortical map reorganization after the training phase, measured by the change in distance between the centroids of the cortical representations of index finger and ring finger, measured in cell units (1 unit = distance between the centers of two neighboring cells on the cortical map), relative to the value on the PRE condition.

(a) Map reorganization for model variation A involving an integrated cortical map receiving input from both tactile and nociceptive channels. (b) Map reorganization for model variation B involving a split map receiving separate input from tactile and nociceptive channels. Positive/negative values correspond to a decrease/increase in distance between the centroids of index and ring finger. Values are expressed as medians, error bars indicate the (25,75)% quantiles, asterisks indicate significant differences to other values (above braces) and to zero (above numbers) (***: p < 0.001, n.s.: not significant, Wilcoxon rank sum test, statistically corrected).

The cortical representation of the amputated finger was preserved after amputation (Figure 4ab, light gray bars on conditions NOPAIN and PAIN are significantly positive).

The cortical representation of the amputated finger during imagined movements of the phantom was significantly more active on the PAIN condition than on the NOPAIN condition (Figure 4b, light gray bars on conditions NOPAIN and PAIN).

Discussion

Our results provide an interpretation of conflicting findings212 concerning the cortical representations in phantom pain by means of two variations of a computational model. The model simulations yielded that the nociceptive channels of the amputated finger showed a significantly stronger central nociceptive activity during the resting phase (no stimulation) on the PAIN condition than on the NOPAIN condition (Figure 4a, dark gray bars). We therefore identify the NOPAIN and PAIN conditions with scenarios where the subject experiences weak and strong phantom pain, respectively. On the basis of this identification, the model predicts that the degree of reorganization in the somatosensory map is stronger in patients with phantom pain than in patients without (or with less) phantom pain. This prediction would be in accordance with the maladaptive reorganization model. Furthermore, however, the model predicts that the representation of the phantom is preserved after amputation, regardless of there being phantom pain or not, and that the activity of the phantom representation during executed phantom movements is stronger in patients with phantom pain than in patients without phantom pain. These latter two predictions would be in accordance with the persistent representation model. Summarizing, in view of our simulations the two explanatory models and their respective experimental findings can be reconciled.

Apart from explaining existing data, the model makes a prediction, which is, to our knowledge, not yet covered by experimental evidence, and which could be easily tested by a questionnaire: Amputation patients suffering from phantom pain should feel fewer phantom sensations that are not painful than patients suffering, or suffering less, from phantom pain (Figure 4a, white bars).

We shall now discuss the individual model assumptions. According to assumption A1, the somatosensory cortex reorganizes itself in response to somatosensory input. This phenomenon has empirically been demonstrated for diverse types of stimulation, including nociception18192021222324. The Kohonen map14 used in our simulations is a rather abstract and idealized model of a biological self-organizing neural network. There are physiologically more elaborate and less idealized models for self-organizing maps encoding the location of receptors and other receptor properties, taking into account also temporal correlations in the receptor input2526. Despite its known limitations, however, the Kohonen map is a numerically efficient and well-established tool to simulate the consequences of amputation and other sorts of sensory deprivation in the sensory cortex272829. In the context of our study, the relevant function of the cortical map is its capability to organize itself according to the topology of the receptor space, which makes the Kohonen map a suitable choice here.

In model variation A, the integrated cortical map only contains where-information, that is, information about the location of the receptor on the skin, and no other receptor properties. In particular, it does not contain what-information, that is, information about the modality of the stimulus, which is here either touch or nociception. Recent findings show that the topological representations of non-noxious tactile stimuli on the skin largely overlap with those of noxious tactile stimuli30, although the representations slightly differ on a smaller scale. Thus, it is not yet resolved, whether or not the modality is encoded in the location of the representation. In any case, the what-information is certainly processed in additional pathways; especially the encoding of painfulness is likely located in the insular-opercular region rather than S13132. The different aspects about a stimulus are presumably at some stage bound together to yield a unified percept, though it is yet not fully understood, and part of the notorious binding problem33, as to how this unification is actually accomplished in the brain. Here, we do not attempt to resolve this problem but rather computationally unify where and what-information (Figures 3–5).

In model variation B there are two modality-specific cortical maps. So, the two maps together encode what-information, in contrast to model variation A. In a biological system, these two modality-specific maps could either be spatially separate, or they could overlap to some extent. In either case, observable differences in cortical activity would strongly depend on the applied measurement method. If cortical reorganization is measured using tactile stimuli only, the model predicts that one observes more reorganization on the PAIN condition than on the NOPAIN condition, in accordance with the results of Flor et al2, supporting the maladaptive reorganization model (Figure 5b, white bars). If, on the other hand, cortical reorganization is measured using nociceptive stimuli only, the amount of reorganization is predicted to appear smaller on the PAIN condition than on the NOPAIN condition (Figure 5b, dark gray bars). This model prediction, which is to our knowledge not yet covered by experimental evidence, could be tested in a manner similar to the procedure used to establish the maladaptive reorganization model2, only with using nociceptive rather than non-nociceptive stimuli.

A third possible model variation, which we have not included in our simulations, is the case of a single unified cortical map that receives input from both tactile and nociceptive channels (as in model variation A), but which organizes itself not only according to the position of receptors (as in A and B), but also according to the modality of the stimulus. There are modeling approaches along this line, with a focus on the optic tectum and the primary visual cortex2526. The resulting cortical structure contains modality-specific regions that are interspersed with each other, similar to the orientation-specific columnar structure of the primary visual cortex3435. In the case of touch and nociception, such regions would, in effect, form two modality-specific sub-maps, a structure that resembles the situation described by our model variation B. While a classical Kohonen map, which has been implemented in our model, is indifferent with respect to temporal correlations between the afferent signals, the emerging structure of the biological somatosensory cortex will probably be influenced by the temporal coherence between nociceptive and tactile events. Therefore it is possible that the topological structure of the real somatotopic map also encodes modality to some degree. The existing experimental evidence is not decisive about this question, and recently measured nociceptive and non-nociceptive representations of fingers in the primary somatosensory cortex turned out to be highly aligned at the resolved scale30.

Figure 2 shows a representative example from the set of simulated map formations. The exact structure of the simulated maps differed to some degree from simulation to simulation, due to the stochastic nature of the initial state of the map and due to the stochastic occurrence of spontaneous events during the training phase. At the end of each simulation, however, the fingers were mapped to a corresponding coherent region on the cortical map, and the spatial relation between the fingers was always represented as well.

Assumption A2 postulates spontaneous neural activity in the sensory system that is increased in those parts affected by deafferentation, and which is implemented in the computational model in the form of discrete neuronal noise (DNN) and spontaneous coherent activity (SCA). The origin of the increase of spontaneous activity is still not completely understood, but there are several lines of explanation. Spontaneous activations of afferent fibers, also known as ectopic discharges, have been found after peripheral nerve lesion and neuropathic degeneration36373839404142434445. In certain cases the increase of spontaneous activity can be traced back to long-term alterations of the nervous system due to previous and recent nociceptive activity, a mechanism that is often referred to as pain memory46474849505152. As for the DNN, increased spontaneous activity of nociceptive fibers at the painful site has been measured in patients suffering from spontaneous pain in connection with CRPS53, diabetic neuropathy38 and phantom pain36 (measured in neuromata at the stump). Moreover, an increased spontaneous activity of nociceptive fibers in rats has been demonstrated to be a consequence of spinal cord injury4345. It has been argued that an increased spontaneous activity of C- and Aδ-fibers is a possible cause for dysesthesia, that is, spontaneous pain4054. As for the SCA, in the model this type of spontaneous activity is implemented to occur after the spinal gate in higher regions of the nociceptive pathway. The SCA is modeled in the form of coherent activity patterns resembling those patterns elicited by actual noxious stimulation. This type of activity may be a possible cause of more detailed painful experiences. About 40% of the patients describe their phantom pain as being close to actual pain experienced in the limb before or during its amputation, and these sensations are often referred to as pain memories515556575859. As for the possible site of the physiological mechanisms leading to pain memories, studies show that not only the cortex but also subcortical regions such as the cerebellum are involved in pain-related associative learning6061. The central role of pain memories in the context of chronic pain and phantom pain has often been emphasized and is a matter of current debate1135662. One might conceptualize the SCA in our model as spontaneous occurrences of detailed pain memories, although it is not necessary to do so.

Assumption A3 postulates that phantom movements and attention to the phantom limb activate neural circuits that have stored sensory experiences of the missing limb. This assumption would be compatible with studies showing a modulation of phantom pain during the perception of the mirrored intact limb at the site of the phantom56. It would also be compatible with studies reporting phantom pain experiences elicited by concentration on the phantom limb58, by watching individuals whose corresponding intact limb is touched9 or by observing pain in others6364. A crucial measurement in the study of Makin et al.12 in support of the persistent representation model, consisted in having the subject perform phantom movements while measuring the activity of the somatosensory cortex. To model an analogous scenario, our simulations involved, subsequent to a training phase where the model is trained with input from random somatosensory stimulations, a probing phase. Given that the subject executes a phantom movement, then according to assumption A3 the movement execution would activate neural circuits that increase the strength of spontaneous coherent activity (SCA) in the somatosensory channels corresponding to the moved (phantom) limb. In our simulations, during the probing phase the SCA strength of the affected somatosensory channels was multiplied by a factor of five, and the resulting activity of the cortical map was calculated.

Assumption A4 postulates that a central gate is regulating the input to the somatosensory cortex. The central gate is probably the most speculative element in our model. A likely location for it might be the thalamus which is known to be a relay station for all afferent connections into the cortex65. There are two more gates in the model, the peripheral gate, and the spinal gate, which are less speculative. The peripheral gates correspond to receptors in the periphery, and the opening and closing of each gate would correspond to the lowering and raising of the gating threshold, leading to peripheral sensitization and desensitization, respectively. The gates are implemented as linear saturated functions (Figure 6). For stimuli of limited strength the linearity of the receptor response is a reasonable approximation66. The cut-off at the channels' maximum excitability has been implemented to take into account that any receptor will go to saturation for strong enough stimuli [ibid.]. Further, the spinal gates are simplified versions of neural mechanisms in the dorsal horn of the spinal cord, which have first been discovered and modeled by Melzack and Wall in the context of their seminal gate control theory151617. Their closing and opening contributes to central sensitization and desensitization, respectively, which is induced and maintained in the body by activity in peripheral afferent fibers4054, by neurogenic inflammation [ibid.], and also by descending supraspinal modulation6768.

Figure 6

Input-output relation of the linear gates implemented in the model.

External modulation affects threshold and gain, the output is cut off at unity.

There are certain limitations of the computational model. The model involves separate modality-specific channels from the periphery up to the central gate before the somatosensory cortex. This is an idealization in so far as there are wide dynamic range (WDR) neurons in the deep dorsal horn (laminae V–VI) that respond to stimuli in both the noxious and non-noxious domain, so that their firing rate encodes the strength of the stimulation but not the modality6970. The functional role of WDR neurons is still controversial, and some researchers hold that the WDR neurons, rather than the nociceptive specific (NS) neurons, are responsible for the subjective perception of the intensity of painful stimuli717273. In any case, we do not expect the results of our numerical simulations to differ qualitatively if WDR neurons were included in the model. Our expectation is based on the fact that the self-organizing map used in the model encodes only where-information (see above), so the map would be reorganized regardless whether or not the altered input exclusively comes from nociceptive-specific channels.

The maladaptive reorganization model and the persistent representation model are each based on an empirically established relationship between a certain sort of physiological fact on one side and a subjective report of pain on the other side. These empirical findings are in so far unproblematic as a correlation is an objective statistical property of given data, while pain reports are an objective (though not necessarily reliable) measurement of a subjective experience. In order to relate our simulation results to the empirical findings, we have to assume a relationship between the simulated physiological state and a subjective experience of pain. We cannot “ask” the model to what extent it is in pain during the simulation. For this study, we have defined the central nociceptive activity as the remaining activity of nociceptive channels after having passed the central gate, and we have taken this value as an estimate for the quantity of subjectively experienced pain. Generally, there will be further modulation by high-level processes, including cognition, visual perception, psychological influences, and learning, which are not covered by the computational model.

The simulations suggest that phantom pain, maladaptive reorganization, and persistent representation may all be caused by the same underlying mechanism. Accordingly, the reorganization of the somatosensory cortex would rather play the role of an epiphenomenon that is correlated with, and therefore acting as a marker for, phantom pain. In our model, however, maladaptive reorganization in S1 would not cause the painful experience itself; instead, the causal driver of phantom pain, as far as our model suggests, is the abnormally enhanced spontaneous activity of deafferented nociceptive channels. As the model lacks a perceptual system, there is, however, room for potential causal influences on phantom pain other than spontaneous activity. A further limitation of the model is the lack of an underlying mechanism for the abnormal enhancement of spontaneous activity after deafferentation, which therefore remains unexplained.

To summarize, simulations of a computational model built upon physiologically plausible assumptions might help to reconcile two apparently contradictory empirical findings and their corresponding conceptual models. In agreement with one of these findings, the computational model predicts that an abnormally increased spontaneous neural activity following amputation induces cortical reorganization that is more pronounced in patients suffering from phantom pain as compared to patients without phantom pain. In agreement with the other one of these findings, the activity of the cortical representation of the missing limb during executed phantom movements is predicted to be stronger in patients suffering from phantom pain as compared to patients without phantom pain.

Methods

The computational model (Figure 1) is programmed in MatLab™ using the SOM Toolbox (Copyright 2000–2005 by E. Alhoniemi, J. Himberg, J. Parhankangas and J. Vesanto), and it is designed as follows. Each finger of a two-dimensional image of a hand is equipped with tactile and nociceptive receptors that are randomly distributed according to a homogenous two-dimensional receptor density ρ, so that there typically are N = Aρ receptors in total, with A being the total surface area of all fingers. On this area, each receptor i has a randomly selected location x = (x, y) and a modality value corresponding to tactile and nociceptive modality, respectively. The receptor is stimulated by a Poisson process of Gaussian-shaped stimulus events of a given duration whose amplitude varies randomly between zero and a given maximum amplitude smaller or equal to unity, occurring randomly at a given average rate of λ events per second (Table 1, parameters “stim dur”, “stim amp”, and “stim rate”, respectively). Each receptor is connected with one and the same cortical map by a neuronal channel that is interrupted by three linear gates. The first gate, the peripheral gate f, corresponds to the receptor's sensitivity towards the stimulus S in such a way that if S exceeds a certain threshold , the receptor generates action potentials with a firing rate that is linearly related with the stimulus strength S(t) by a gain factor g1 ≥ 0. The output of the gate is cut off at unity, which guarantees that the firing rate remains within the range [0,1]. To reach the somatosensory cortex, the signal generated by the receptor has to pass two more linear gates, the spinal gate f2 and the central gate f3, each endowed with their own threshold and gain. The general formula for all three gates f reads where x is the respective input to the gate (Figure 6). For each gate, the default gain in our model simulations has been chosen in dependence of the default threshold of the gate, so as to approximately preserve the mean strength of the input signal corresponding to zero threshold. Thus, for the default threshold of θ = 0.1 (i = 1,2,3), the gain was set to (Table 1).

After passing the peripheral gate f1, discrete neuronal noise (DNN) is added to the signal in form of a Poisson point process N(t) with an average rate of λ point-like spikes per second, with the spike amplitude varying randomly between zero and a certain maximum value (below unity) (Table 1, parameters “DNN rate” and “DNN amp”, respectively). Next, the signal passes the spinal gate f2, after which spontaneous coherent activity (SCA) is added to the signal in form of a Poisson process M(t) with a given average rate λ. In contrast to the discrete neuronal noise N(t), the SCA are not spikes but rather bursts, which are modeled as Gaussian packets of a given duration and a given amplitude (Table 1, parameters “SCA dur” and “SCA amp”). The signal finally has to pass the central gate f3, so that after the central gate the firing rate yields The occurrence times of external stimulation S(t), neuronal noise N(t), and spontaneous coherent activity M(t), are modeled as Poisson processes with rates λ, λ, λ, respectively, each corresponding to the average number of events per second, so with Δt being the size of a simulated timestep, the probability for an event to occur within one time step would read p = λΔt. The simulations were run with a temporal resolution of Δt = 0.1 simulated seconds.

After the central gate, the signal enters the cortical map, which is modeled as a self-organizing Kohonen map14. It is a rectangular grid of neurons receiving input from the channels, and the synaptic weights of this neural network are adapted by a certain update rule taking into account the topological structure of the map as well as the receptor location. Initially, the input weights of the neurons that constitute the map are set to random values. In our simulations, the map was trained in a batch to increase numerical efficiency. In batch training, all channel activations passing the central gate during training time, are presented to the map at once in sequential iterations n = 1,…,N. Each input (channel activation) is of the form x = (x,y), where x and y are the coordinates of the location of the receptor that corresponds to the activated channel. Let xs,s = 1,…,S, be the list of all inputs occurring during training time, with S being the total number of inputs, then the weight update rule is given by where (n + 1) is the weight vector of the i-th map unit at iteration step n + 1, where is the neighborhood function between the map unit i and the best matching unit (BMU) c corresponding to the input x. The neighborhood function is given by where is the distance function for two neurons i and j located at positions and on the grid, where σ(n) is the radius of the neighborhood. The BMU c corresponding to input x is defined by with being the Euclidean norm. The map was trained in a batch with two phases: In the “rough” phase, the iteration length is set to N = 50 and the radius σ decreases linearly from 20 to 5; in the “fine” phase the iteration length is set to N = 20 and the radius σ decreases linearly from 5 to 1. During training, the cortical map forms a stable somatotopic (neighborhood-preserving) representation of the area over which the receptors are distributed. Hence, if receptor i at location x = (x,y) is stimulated, and the signal passes all three gates, then a certain neuron c in the cortical map, the BMU, is activated.

However, not every neuron is a best matching unit to some receptor, so certain units would remain “blank”. If neuron c in the cortical map is activated by a signal of strength r in a tactile or nociceptive channel, and the synaptic weights of c encode the location (x,y) on the hand, then the representation of this activation is respectively a tactile or nociceptive stimulus of strength r at the location (x,y). These representations might correspond to tactile or painful percepts induced by stimuli or spontaneous neural activity, if areas of the brain associated with conscious perception, which are not included in our model, would further process the signals.

The numerical simulation is divided into a training phase, a probing phase, and a resting phase, on each of three conditions denoted as PRE, NOPAIN, and PAIN, corresponding to different parameter values representing different scenarios (Table 1).

The PRE condition corresponds to a healthy subject with no amputated fingers. The NOPAIN condition corresponds to the same subject after amputation of the middle finger, so any stimulation of the channels originating in the amputated finger is suppressed. Also, the thresholds of the spinal and central gates of the affected channels are lowered to 1/4 of the PRE value, which is conceptualized as being a consequence of the deafferentation. The PAIN condition also corresponds to a post-amputation scenario; this time, however, there is a strong enhancement of spontaneous activity in the channels affected by amputation: the SCA amplitude and frequency are increased by a factor of five as compared to the PRE condition. Also, the central gates of the affected channels raised to 150% of the respective PRE values, which is conceptualized as being a consequence of the strong SCA, while the spinal gates are lowered to 1/4 of their PRE value. While the simulation of the PRE condition started with an initially randomized cortical map, the simulations of the NOPAIN and PAIN condition were each based on the same cortical map that resulted from the simulation of the PRE condition.

In the training phase on the PRE condition, the cortical map was set to a random initial state, all receptors on the hand were stimulated with randomly occurring Gaussian shaped events of constant width (Table 1, parameter stim dur) and random amplitude between zero and a certain maximum value (Table 1, parameter stim amp), the channels were processed, and the map was trained with the channel output. In the training phase on the NOPAIN and PAIN conditions, the cortical map was set to the state resulting from the training phase on the PRE condition, all receptors except those located on the middle finger were stimulated, the channels were processed, and the map was trained with the channel output. After the training phase, the cortical map showed a somatotopic organization, with each finger corresponding to a region on the cortical map with preserved neighborhood relations (Figure 2).

The training phase on each condition is followed by a probing phase, where the parameters are set to the same values as in the training phase, except that the SCA rate and amplitude further multiplied by a factor of five. The multiplication is conceptualized as resulting from a voluntary movement of the phantom finger during the probing phase. During the probing phase, the map is not trained but kept constant.

The final resting phase on each condition corresponds to a phase without any stimulation, with the cortical map being in the same state as in the probing phase, and with the model parameters being set to the same values as in the training phase. The cortical map is not trained during this phase.

The simulated durations of the phases are: training phase: 60 s; probing phase: 240 s; resting phase: 300 s. The duration of the training phase was chosen to be long enough for the cortical map to become stable, which has been tested in separate simulations. The durations of the probing and resting phases were chosen such that there was enough accumulated activity in the map to yield sufficiently stable results, in particular for the (rather weakly active) nociceptive channels.

The cortical map in Figure 2 was obtained by activating one by one all channels originating in a certain finger, and coloring those cells in the cortical map that are the best matching units corresponding to the location of the receptor that belongs to the channel.

The activation map on the left side of Figure 3 was obtained by accumulating the activity (in terms of normalized firing rate) of each cell in the cortical map during the probing phase. The representations of the cells depicted on the right side of the figure were obtained by reading out the positions encoded by each of the activated cells in the cortical map, and by giving them the same color as the cell.

Figure 4 depicts values resulting from thirty simulations of model variation A. The corresponding values for model variation B are only statistically different, so they are omitted here. On panel (a), the values for the central nociceptive phantom activation were obtained by summing the activity of all channels that originate in the missing finger and accumulate them over the resting phase. On panel (b), the normalized firing rate of all units of the map was summed and accumulated over the probing phase, to obtain a numerical measure of cortical map activity.

In Figure 5, the degree of cortical map reorganization on the NOPAIN and PAIN conditions with respect to the PRE condition was measured by r = (d − d), where d and d are the distances between the centroids of the representations of the index finger and the ring finger on the PRE and POST condition, respectively, and where POST is either NOPAIN or PAIN. Panel (a) shows the map reorganization for model variation A involving an integrated cortical map that receives input from both tactile and nociceptive channels. Panel (b) shows to the map reorganization for model variation B involving two separate modality-specific cortical maps that receive input from tactile and nociceptive channels, respectively.

The computational model contains stochastic elements, such as the randomly initialized cortical map, the discrete neuronal noise (DNN), and the spontaneous coherent activity (SCA). Thus, to obtain reliable predictions we performed thirty simulations for each model variation A and B, corresponding to integrated and split cortical map, respectively. The resulting values of interest, which are depicted in Figures 4 and 5, may not be normally distributed, in particular since some of them are bounded from below by zero. Thus, we have taken the median as the main value and the (25,75)% quantiles as error bars in Figures 4 and 5. Furthermore, instead of t-tests we used Wilcoxon rank sum tests (equivalent to Mann-Whitney U-tests) to check for statistical significance. The resulting p-values have been statistically corrected using the Bonferroni method. Although the Bonferroni method is very conservative, all significances but one (the reorganization of the nociceptive map on the PAIN condition) turned out to be very robust, yielding p-values far below the limit of p < 0.001. We have run significance tests for any combination of differences, including the difference to zero baseline, without loosing high significance (p < 0.001). For the sake of clarity, though, in the Figures we have only indicated the significant differences (or the lack thereof) necessary to support the central statements made in the Results section.

Author Contributions

K.J.B. and H.W. programmed the computational model. K.J.B. wrote the manuscript. M.d.L., H.W., C.P. and T.W. commented on and reviewed the manuscript. K.J.B. and M.d.L. prepared the figures. All authors contributed to the conceptual foundations of the model.