Participates

Study participate criteria: 1. Meet the diagnostic criteria for disorders caused by DSM-5 psychoactive substances; 2. Patients within six months of the withdrawal period; 3. Junior high school education and above; 4. Years of Age 18–41; 5. Voluntary participation in this study and sign the informed consent form. Exclusion criteria: 1. Severe cognitive dysfunction, unable to cooperate with the completion of project-related assessment and testing; 2. Patients with severe physical diseases; 3. Severe psychotic symptoms; 4. Have other mental activities Substance abuse (except nicotine). The study was conducted in accordance with the declaration of Helsinki and was approved by the Ethics Committee of Shanghai University (Approval No. ECSHU2020–071).

Drug users description

According to the demographic scale made in the early stage, the most important items of 30 people (All male), such as drug use type, drug history, average drug dosage and drug use frequency, were accurately inquired and counted.

Methamphetamine easily causes intense excitement, which is difficult to eradicate after addiction.

Heroin is a psychoanaesthetic drug. Once a person becomes addicted, their physiological reaction is intense, and they have a compulsion to seek medication.

“Mixed drug abusers” refers to drug users who attempt to mix two or three drugs at a time. The main types of drugs used are: LSD, Flunitrazepam, N2O, Pethidine hydrochloride, MDMA, Cannabis, ketamine, etc. Table 1 is the personal information of the selected 30 subjects.

Table 1 Personal information of drug addictsNIRS technology equipment introduction

This paper uses a high-density NIRS device (NIRSIT; OBELAB, Seoul, Korea), and its specific hardware parameters and wearing methods are as follows: the specific hardware parameters are shown in Table 2. Figure 1 NIRSIT wearing method in the experiment.

Table 2 NIRSIT specific hardware parametersFig. 1

NIRSIT wearing method in the experiment

NIRS channel and functional area division

The four advanced functional areas detected by the forehead fNIRS device are divided into the dorsolateral prefrontal cortex, the ventrolateral prefrontal cortex, frontopolar prefrontal cortex, and the orbital frontal cortex. The specific channel distribution: the right dorsolateral prefrontal cortex is 1,2,3,5,6,11,17,18 channels. There were 19, 20, 33, 34, 35, 38, 39 and 43 channels in the left dorsolateral prefrontal lobe. There are 4,9,10,40,44,45 channels in ventrolateral prefrontal cortex of left and right hemispheres. There are 14, 15, 16, 29, 30, 31, 32, 46, 47, 48 channels in the left and right orbital frontal cortex. Frontopolar prefrontal cortex is 7,8,12,13,21,22,23,24,25,26,27,28,36,37,41,42 channels. Fig. 2 NIRS channel and functional area division.

Fig. 2

NIRS channel and functional area division

Near-infrared imaging theory

When light passes through a uniform, non-scattering medium, only the absorption effect of the medium on the photons is considered. According to the Beer–Lambert law [34,35,36,37], the attenuation of light intensity is expressed as follows:

$$ OD=\log \frac=-\upvarepsilon \left(\lambda \right) cd\kern0.1em \log \kern0.3em \mathrm $$

where I0 is the incident light intensity, I is the incident light intensity, ε(λ) is the extinction coefficient of the substance at a wavelength of λ, determined by the absorbing medium and the wavelength of the light, c is the medium concentration, and d is the thickness of the medium. The absorption coefficient μa is defined as follows:

$$ _a=\varepsilon \left(\lambda \right)c $$

The total absorption coefficient of the medium can be expressed as a linear superposition of the absorption coefficients of each medium:

$$ _a\left(\lambda \right)=\sum \limits_i^N_i\left(\lambda \right)_i $$

The optical density is the product of the medium thickness d and the total absorption coefficient μa [38, 39]:

$$ OD=\mathit\frac=d\sum \limits_i^N_i\left(\lambda \right)_i $$

The actual biological tissue is very complex and is a strong light scatterer, and light undergoes multiple scatterings during its output before it can be detected. The attenuation of light in tissues includes absorption and scattering, and to capture the effect of scattering on light loss, Deply et al. proposed a modified Beer–Lambert law [40], expressed as follows:

$$ OD=\log \frac= DPF\left(\lambda \right)\cdot dc\upvarepsilon +G $$

where G denotes the light loss due to scattering and other boundary losses. DPF is the differential path factor, whose value is the ratio between the actual optical path length traveled by light in the tissue and d. The DPF values in different tissues can be obtained from the literature [41].

Measurement of changes in hemodynamic parameters

To detect changes in hemodynamic parameters using NIRS, a reference state is usually selected in near-infrared measurements to detect relative changes in the concentration of absorbing chromophores Δc, based on the modified Beer–Lambert law [42, 43]:

$$ \varDelta OD=\mathit\frac= DPF\left(\lambda \right)\cdotp d\varDelta c\varepsilon $$

When detecting relative changes in HbO2 and HHb concentrations:

$$ \varDelta O^=\left(_^\varDelta \left[ Hb_2\right]+_^\varDelta \left[ HHb\right]\right) DPF\left(\lambda \right)\cdotp d $$

where Δ[HbO2] and Δ[HHb] are the variations in HbO2and HHb concentrations, the selected incident near-infrared wavelengths λ1 and λ2, brought into the above equation have [44, 45].

$$ \varDelta O^=\left(_^\varDelta \left[ Hb_2\right]+_^\varDelta \left[ HHb\right]\right) DPF\left(_1\right)\cdotp d $$

$$ \varDelta O^=\left(_^\varDelta \left[ Hb_2\right]+_^\varDelta \left[ HHb\right]\right) DPF\left(_2\right)\cdotp d $$

When the value of DPF is known, solving the above equation for the system of equations yields the change in HbO2 concentration Δ[HbO2] and the change in HHb concentration Δ[HHb], expressed as follows:

$$ \varDelta \left[ Hb_2\right]=\frac^\frac^}_1\right)}-_^\frac^}_2\right)}}_^__2}^-_^__2}^\right)} $$

$$ \varDelta \left[ HHb\right]=\frac_2}^\frac^}_1\right)}-__2}^\frac^}_2\right)}}_^__2}^-_^__2}^\right)} $$

Experiment design and data acquisition

We used E-prime software package (Psychology Software Tools, Pittsburgh, PA) to write the experimental paradigm, with each map numbered. A complete experimental paradigm consists of the following three stages.

The first stage of the experimental paradigm, 10 min in total, during which the subjects need to close their eyes for 5 min and then open their eyes for 5 min.

The second stage, it lasts 6 min and is divided into drug maps and neutral maps. Among them, each block lasts 10 s. There are a total of 16 maps, and the display time of each map is 0.6 s. At the beginning, the first four maps are displayed randomly, during which there are two drug maps. After displaying the first four maps, the remaining 12 neutral maps are displayed randomly. After a block ends, there will be a 4-s interval map with a white background and a black cross. Figure 3, the examples of drug abuse-related maps used in the experimental paradigm.

Fig. 3

Examples of drug abuse-related maps in the experimental paradigm

The third stage, it lasts a total of 4.6 min, during which the maps are all neutral, with each block lasting 10 s. There are 16 maps in total, with a display speed of 0.6 s. There will be a 4-s interval between each block. Figure 4, the examples of neutral maps used in the experimental paradigm. Figure 5, The whole process of experimental paradigm.

Fig. 4

Examples of neutral maps in the experimental paradigm

Fig. 5

The whole process of experimental paradigm

Linear discriminant analysis algorithm principle

Linear discriminant analysis (LDA) is a supervised pattern recognition method. The LDA classifier reduces the dimensions of the data, reduces complex features into low-dimensional features through projection, and searches for specific classification surfaces to maximize the discrimination between the two types of task classification in order to realize feature classification [46].

LDA is the most traditional linear classifier; there are k linear functions for a k-classification problem.

When yk > yj for all j, then x belongs to class k. When k = 2, it becomes a binary problem.

Support vector machine algorithm principle

Support vector machine (SVM): For nonlinear feature input, we design mapping rules to map nonlinear features to high-dimensional space and ensure that the features are distributed linearly in high-dimensional space. Then, we can easily construct the optimal classification hyperplane in high-dimensional space based on the structural risk minimization criterion so that this optimal classification hyperplane can correctly separate as many of the two types of samples on the one hand and maximize the classification interval between the two types on the other [47].

The sample set is designated as (xi, yi), i = 1, 2, ⋯, l, xi ∈ Rd, where yi ∈  is the category number. The linear separable and linear non-separable cases are explored jointly, and the relaxation variable is introduced ξi ≥ 0, where ξi = 0 represents linear separability, and ξi > 0 represents nonlinear separability. If the classification surface equation is w · x + b = 0 (w is the weight vector, and b is the offset), the classification interval is equal to \( \frac \). Maximizing the classification interval is equivalent to minimizing ‖w‖ or (‖w‖2). To make the classification surface classify all samples correctly as much as possible, it is necessary to meet the following constraints:

$$ _i\left[\left(w\cdotp _i\right)+b\right]-1+_i\ge 0,i=1,2,\cdots, l $$

Therefore, if the constraint formula is satisfied, the (‖w‖2) minimum classification surface can be the optimal classification surface.

Machine learning algorithm classification

In this paper, we propose and design a neural network algorithm model that extracts the hidden features in the 0.625-s fNIRS data after drug map stimulation by convolution to define the type of drugs of abuse. In this paper, a convolutional neural network (CNN) model is designed. It is used to implement the classification of people who abuse different type of addictive drugs. The Structure of CNN used in this paper is shown in Fig. 6.

Fig. 6

Structure of convolutional neural network

Data preprocessing

Butterworth filter, the frequency response curve in the passband has a relatively flat and undulating character, and gradually decreases to zero at the edge of the stopband. In this paper, a Butterworth filter based on infinite impulse response is selected to band-pass filter the acquired Near infrared signals for the purpose of physiological artifact removal. The expression for the n-order Butterworth filter is as follows:

$$ \left(\boldsymbol\right)\right|}^}=\frac}+}}_}}\right)}^\boldsymbol}}=\frac}+}^}}}_}}\right)}^\boldsymbol}} $$

Where n is the order, fc is the cutoff frequency, and fp is the passband edge frequency.

In this paper, n is 6 and the frequency band range is 0.01 Hz to 3 Hz. This band range can remove the interference of heartbeat respiration and slow drift to the raw data, and also can maximize the preservation of hemodynamic characteristics.

The formula for each structure in the CNN network is as follows:

Convolutional layer

Convolutional layers are the core of convolutional neural networks [48]. The calculation form is as follows:

$$ _j^l=f\left(\sum \limits__j}_i^\cdotp _^l+_j^l\right) $$

\( _j^l \) is the jth feature of the layer l.\( _^l \) is the jth feature of the layer l and the ith feature of the layer l − 1. \( _j^l \) is a bias parameter, f(•) is the activation function.

Pooling layer

The pooling layer sub-samples the input features according to specific rules in order to make the network robust to small changes in previously learned features [49]. The calculation form is as follows:

$$ _j^l=f\left(_1^l\backslash down\left(_j^\right)+_j^l\right) $$

\( _j^l \) is the jth feature of the layer l.\( _1^l \) is the Subsampling coefficient.\( _j^l \) is the bias parameter, down(•) is a sub-sampling function, f(•) is the activation function.

Normalization of data

Batch standardization layer: in training convolutional neural network, the input data are usually whitened, which can speed up the training speed.

Set the data value input =  of the input data block, the parameters to be learned are γ and β, first calculate the average value of each data block:

$$ _B=\frac\sum \limits_^m_i $$

Calculate the variance of each data block:

$$ _B^2=\frac\sum \limits_^m_i-_B\right)}^2 $$

Normalize each set of data:

Using the parameters that need to be learned in the network and linear transformation:

$$ _i=\gamma }_i+\beta \equiv B_\left(_i\right) $$

Activation function

In this paper, the activation function uses a modified linear unit (ReLU), and the corresponding calculation formula is as follows:

$$ \mathit}L\mathrm=\left\x,x\ge 0\\ x,x<0\end\right.=\mathit\left(0,x\right). $$

The results show that the derivation of the activation function is simple, and the output of some neurons is 0, which realizes the sparsity of the network, reduces the interdependence of parameters [50].

Full connection layer

Each feature must be converted to one-dimensional before it can be used as the input of the fully connected layer [51]. The calculation as follows:

$$ _(x)=\uptheta \left(^}x+b\right) $$

hw, b(x) is the output value of the neuron. x is the input feature vector of the neuron. w is the weight. b is the bias parameter. θ(∙) is the activation function; The first fully connected layer in this paper uses the ReLU activation function.

Softmax layer

In CNN, if the final output result is single-label multi-classification, the softmax function is usually used to normalize and map to the probability value, and the Softmax calculation formula is as follows:

$$ _i=\mathrm\left(_i\right)=\frac_i\right)}_c\right)} $$

oi is the value of the output neuron corresponding to the Ith category.