I am trying to generate a set of equations that describe a change in density for each variable in the form of: $\frac{dX_i}{dt}=X_i(r_i-s_i*X_i+\sum_{j=1,j\neq i}^na_{ij}X_j)$ where n is the number of variables/species. The code I have generates random values for $a_{ij}$ and produces a list of right hand sides of the equations I need. Next I want to solve them where the right hand sides equal to 0

nSpecies=5;    
S = -0.5;
M = Table[0, {i, 1, nSpecies}, {j, 1, nSpecies}];
Do[Do[
  If[i == j, M[[i, j]] = S ,
   If[RandomReal[] < 0.6,
    M[[i, j]] = -RandomVariate[HalfNormalDistribution[1]];
    M[[j, i]] = -RandomVariate[HalfNormalDistribution[1]];
    ,
    M[[i, j]] = RandomVariate[HalfNormalDistribution[1]];
    M[[j, i]] = RandomVariate[HalfNormalDistribution[1]];
    ]
   ],
  {j, 1, i}], {i, 1, nSpecies}]
(M) // MatrixForm;
    
variables = Take[Alphabet[], nSpecies];
IM = variables.M;
eqs = Rationalize[(IM + 1.2)*variables];
Solve[AllTrue[eqs, # == 0 &], variables, Reals]

However, this produces the error "Solve:a is not a valid variable". I have not been able to work out how to fix this. I also have a feeling I might be approaching this incorrectly. any help is greatly appreciated