1
The regression $y = a + bx + \epsilon$ is not symmetric for replacing y and x, as you already suggested via the link (symmetry of linear regression).
However, it is not trivial whether the expectation values, $\hat{a}$ and $\hat{b}$, are going to be different from a symmetric $E(\hat{a}_{xy})=E(\hat{a}_{yx})=0$ and $E(\hat{b}_{xy})=1/E(\hat{b}_{yx})=1$. Possibly the difference in the regression $y \sim x$ vs $x \sim y$ cancel over all possibilities (we will see this is true for $a$). For this we need to "see" all the possibilities and how they average/cancel.
2
We will show that the intercept is symmetric
$$E(\hat{a})=0$$
however the slope is not, and will be smaller than 1
$$E(\hat{b})<1$$
We show this by comparing the non-symmetric regression with a symmetric regression.
3
As symmetric regression we use the line in between the two regression lines lines defined by the y~x and x~y.
In algebraic terms:
$$\hat{b}_{sym} = \frac{\hat{b}_{xy} + \hat{b}_{yx}^{-1}}{2}$$
and
$$\hat{a}_{sym} = \frac{\hat{a}_{xy} - \hat{a}_{yx}\hat{b}_{yx}}{2}$$
4
Note the following about the symmetric slope.
The slope will be symmetric around the angle of 45 degrees
ie, a probability for an angle with $45+\alpha$ is equal to a
probabilty for an angle with $45-\alpha$, also the expectation value
for its angle coefficient will be 1, since:
$\hat{b}_{sym}(x,y)=\hat{b}_{sym}(y,x)$, thus
$E(\hat{b}_{sym}(x,y))=E(\hat{b}_{sym}(y,x))$
(note that we do not have the same for the slope of the asymmetric case $\hat{b}_{xy}(z_1,z_2)=\hat{b}_{yx}(z_2,z_1)$ )
- yet we must also have the symmetry $E(\hat{b}_{sym}(x,y))=1/E(\hat{b}_{sym}(y,x))$, the line that is described by the expectation of the coefficient, if we change the labels xy than we expect a mirror image which has a different b coefficient (the inverse).
- The slope $\hat{b}_{xy}$ will be smaller due to regression to the mean.
Therefore... the symetric slope will be $E(\hat{b}_{sym})=1$, however then since always $\hat{b}_{xy} \leq \hat{b}_{sym}$ we will have $E(\hat{b}_{xy}) \leq E(\hat{b}_{sym})=1$. And the equality is only true when there are no error terms.
See the following image to recall the regression lines $xy$, $yx$ (red curves) and the symmetric one in between (green curve), this image is from the example code below:

If we plot the angle of the symmetric line (green) and the regression line for the y~x model (the red one that has smaller angel with the x-axis) then we see that the angle of the symmetric line distributes evenly around 45 degrees, however the angle of the non-symmetric line $y~x$, is always lower:

5
What is left is to show that $E(\hat{a})=0$. I do not have a direct proof for this. However, it should be sufficient to notice that the regression to the mean (which is zero) does not favor either larger or smaller $a$. Note that there is a point-symmetry, $(x,y) \rightarrow (-x,-y)$, for which the angle $b$ is invariant, but not the intercept $\hat{a}$, which means that it, $a$, should be zero.
You obtain values of $\hat{a}$ different form zero, but you should take into account it's sample variance. (The sample error for $\hat{a}$ is different from zero)
6
Below is adjustment to your code to examine the above described effects
n_pop <- 100
n_sample <- 20
draws <- 1000
M <- matrix(c(1:(2*draws)), nrow=draws)
K <- matrix(c(1:(2*draws)), nrow=draws)
ma <- matrix(c(1:(2*draws)), nrow=draws)
meanx <- matrix(c(1:(1*draws)), nrow=draws)
for (j in 1:draws) {
x <- sort(sample(-n_pop:n_pop,n_sample,replace=FALSE))
y <- sort(sample(-n_pop:n_pop,n_sample,replace=FALSE))
xy_model <- lm(y~x)
yx_model <- lm(x~y)
M[j,] <- coef(xy_model)
K[j,] <- coef(yx_model)
# xy and yx model coefficients last loop
b1 <- coef(xy_model)[2]
b2 <- 1/coef(yx_model)[2]
a1 <- coef(xy_model)[1]
a2 <- -coef(yx_model)[1]*coef(yx_model)[2]
ma[j,1] <- 0.5*sum(a1,a2)
ma[j,2] <- 0.5*sum(b1,b2)
meanx[j,1] <- mean(x)
}
# means
mean(M[,1])
mean(M[,2])
mean(ma[,1])
mean(ma[,2])
# plotting last for-loop x-y as example of regresion fits
plot(x,y)
xp <- -n_pop:n_pop
lines(xp,a1+b1*xp,col=2)
lines(xp,a2+b2*xp,col=2)
lines(xp,0.5*a2+0.5*a1+0.5*b2*xp+0.5*b1*xp,col=3)
# plotting coefficients symmetric versus xy
plot(180/pi*atan(ma[,2]),180/pi*atan(M[,2]), xlab="symmetric fit angle", ylab="xy-fit angle")
lines(c(0,90),c(0,90),col=2)
# plotting mean of x vs beta
plot(meanx,ma[,2]-M[,2], xlab="mu_x",ylab="beta_sym - beta_xy",pch=21,
col=adjustcolor("black",alpha.f=0.1),
bg=adjustcolor("black",alpha.f=0.1),log="y")