For the 1D, 1-channel case, you may be interested in a related question and answer here.
In the 2D case*, if we flatten the input to the convolution $x \in \mathbb{R}^{C\times H\times W}$ into a vector $x' \in \mathbb{R}^{CHW}$ in the usual manner (such that $x'_{iHW+jW+k} = x_{i,j,k}$), and we have a convolutional kernel $K \in \mathbb{R}^{D\times C\times P\times Q}$ ($D$ is the out dimension and each filter has receptive field $P$ by $Q$) then we can define a weight matrix $M \in \mathbb{R}^{DH'W' \times CHW}$ ($H'=H-P+1, W'=W-Q+1$)such that the flattened version of $y = \text{conv}(K,x)$ can be written as $y' = Mx'$ as follows:
$$
M_{s,t} = \begin{cases} K_{s,i,v-j,w-k} &\text{if } 0 \leq v-j < P \text{ and } 0 \leq w-k < Q\\ 0 &\text{otherwise } \end{cases}
$$
Where $i,j,k$ are defined by $t = iHW+jW+k$ and $j<H$, $k<W$. and $u,v,w$ are defined by $s = uH'W'+vW'+w$ and $v < H'$ and $w < W'$.
You can see in each row of $M$, corresponding to a single entry in the output feature map, the only nonzero entries of that row are in the columns corresponding to the appropriate input receptive field.
*well I'm not masochistic enough to deal with strides, dilation, padding, separable filters, etc in this answer.
However, what are neurons in this case?
A neuron in a convolutional network (although I think it's usually not useful to think in terms of neurons), is a single entry in a feature (which is a vector) in a feature map (which is a 2D grid of features -- a 3D tensor).
Ok I agree the indexing notation is rather dense, here I'll write out an explicit example:
Our input $x$ is 1 by 3 by 3:
[
[
[1 2 3]
[4 5 6]
[7 8 9]
]
]
Each single value such as "1" or "5" here is a neuron.
We flatten this into the vector $x'$:
[1 2 3|4 5 6|7 8 9]
(to keep things sane, here and later, i use | to delimit every 3 elements, so you can see how they map onto the 1x3x3 input)
Meanwhile our kernel $K$ is 2 by 1 by 2 by 2:
[
[
[
[a b]
[c d]
]
]
[
[
[e f]
[g h]
]
]
]
We arrange this into the matrix $M$, which is 8 by 9:
[
[a b 0|c d 0|0 0 0]
[0 a b|0 c d|0 0 0]
[0 0 0|a b 0|c d 0]
[0 0 0|0 a b|0 c d]
[e f 0|g h 0|0 0 0]
[0 e f|0 g h|0 0 0]
[0 0 0|e f 0|g h 0]
[0 0 0|0 e f|0 g h]
]
Then $Mx' = y'$ computes
[1a+2b+4c+5d, 2a+3b+5c+6d, 4a+5b+7c+8d, 5a+6b+8c+9d, 1e+2f+4g+5h, 2e+3f+5g+6h, 4e+5f+7g+8h, 5e+6f+8g+9h]
Again as before, each scalar value here such as "1a+2b+4c+5d" is a single neuron.
We reshape this to 2 by 2 by 2 to recover $y$:
[
[
[1a+2b+4c+5d, 2a+3b+5c+6d]
[4a+5b+7c+8d, 5a+6b+8c+9d]
]
[
[1e+2f+4g+5h, 2e+3f+5g+6h]
[4e+5f+7g+8h, 5e+6f+8g+9h]
]
]
And you can see by inspection this is what we'd get by sliding the filter $K$ over the original $x$.