Well, for me this is so obvious that there is no need to answer it. If evolution is not random what else can it be? What is the opposite to random? There are not a lot of alternatives. I don't answer this in the article because I want the reader to answer it for him/herself.

In case someone is interested, I published an update on ResearchGate. Do you know a journal that would be interested in publishing this kind of stuff?

After a bit of reflection, I see what's still missing from my story: everything I wrote is correct, but only for a process that does not terminate. That's where the nuance lies. The problem is now to find the probability distribution for the first time the process achieves a given state (which I haven't calculated yet). At the moment I don't have time to go into depth about this (and I'd have to figure this out for myself as well), but I suggest you take a look at FirstPassageTimeDistribution in Mathematica.

For example, if you have a single letter drawn from A, T, C and G that mutates with probability m, then Mathematica gives you as the probability for achieving the state A for the first time in v steps:

PDF[
  FirstPassageTimeDistribution[
     DiscreteMarkovProcess[{1, 1, 1, 1}/4, transitionMatrix[m]], 
     1
  ], v] = (9 (1 - m/3)^v m) / (4 (-3 + m)^2) 

(for v > 1). I suggest you play around with that.

Thanks for the clarification, that helps a lot. I can now say safely that the problem you're interested in (without natural selection) is definitely a Markov chain (or actually a set of Markov chains). So for this problem, I highly recommend reading up on those, there is a lot of knowledge about them that you can find and will help you.

So if I understand correctly, in the problem without natural selection, every generation each letter (I'm going to stick with letters for now) has a probability to change. Each change (if it occurs) draws a new value from a set of letters of size a completely randomly. From my point of view, this means that the probability distributions of each of the letters are completely independent of each other: all you need to know is the probability P(v) that one of the letters has the right value at a certain generation v. Then you just raise that probability to the power n of letters that you have. Do you agree on that?

This probability P(v) is straightforward to calculate from the theory of Markov chains. I'll try to do the translation as best as I can. Suppose we have as our letters A T G and C and we are going to follow what happens to a particular nucleotide as it mutates. At any given generation v, we will denote the probabilities that the nucleotide is a particular flavor by a vector:

{PA(v), PT(v), PC(v), PG(v)},

where:

PA(v) + PT(v) + PC(v) + PG(v) == 1

(because it's definitely one of the four). So e.g., PA(v) is the probability that the nucleotide has value A at generation v. The vector

{PA(v), PT(v), PC(v), PG(v)}

is called the state vector of the Markov chain at step v. Next, we want to do something with this vector: we want to advance it one step of time. This we do with a matrix that contains the transition probabilities to go from one state to another. For our 4-letter alphabet, this matrix is:

      transitionMatrix[m_] = 
        { 
          {1 - m, m/3, m/3, m/3},
          {m/3, 1 - m, m/3, m/3},
          {m/3, m/3, 1 - m, m/3},
          {m/3, m/3, m/3, 1 - m}
        }

Take the first row of this matrix: these are the transition probabilities if the nucleotide is an A right now. The 1 - m entry (which is on the diagonal of the matrix) is the probability that the nucleotide doesn't change: if it was an A, it stays a A. The off-diagonal m/3 entries are the probabilities that it changes into T (2nd entry in the row: m/3), C (3rd entry in the row: also m/3) or G (4th entry in the row: again m/3). The remaining rows are the transition probabilities for when the nucleotide is currently T, C or G respectively. Note that every row sums to 1.

Now here comes the crucial bit: to get the state

{PA(v+1), PT(v+1), PC(v+1), PG(v+1)}

for the next generation, you perform a left-side matrix multiplication:

 {PA(v+1), PT(v+1), PC(v+1), PG(v+1)} =  
   {PA(v+1), PT(v+1), PC(v+1), PG(v+1)} . transitionMatrix[m] 

Judging from what you wrote down, you already understand the basic idea here, though make sure you really understand this step.

Bear with me for a bit longer, we're nearly there. You assumed the letters at the first generation (v=1) are random, so {PA(1), PT(1), PC(1), PG(1)} = {1/4, 1/4, 1/4, 1/4} Note that this is a left-eigenvector of transitionMatrix[m] (with eigenvalue 1):

{1/4, 1/4, 1/4, 1/4}.transitionMatrix[m] == {1/4, 1/4, 1/4, 1/4 }

This means that the probability vector doesn't change: at every step, the probabilities for A, T, C and G remain 1/4, no matter the value of m!

If you now want the state vector at a given v for any given starting state, you use DiscreteMarkovProcess in Mathematica. For example, try:

FullSimplify[
 PDF[
  DiscreteMarkovProcess[{1,0,0,0}, transitionMatrix[m]][v], 
  k]
]

Here, {1,0,0,0} is the starting state (the nucleotide starts as A with probability 1). v is the generation again, and k is the index of the letters: k = 1 represents A, k = 2 represents T, k = 3 represents C, and k = 4 represents G. If you evaluate this code, you'll see that you get a probability distribution that exponentially decays to the equilibrium state vector {1/4, 1/4, 1/4, 1/4} with increasing v.

So there you have it. I think this should answer your problem, though if I misinterpreted part of the question you should let me know. At the very least this should put you on the right track, I hope.

I get the impression that we don't really understand each other here. Could you repeat your basic problem again, but as clearly as possible? I'm not quite sure I understand what you're actually interested in and what you're doing.

The problem you describe sounds like two Markov chains (or discrete-time Markov processes). Each number (call them x and y) occupy a certain state and each step of the process they have a probability m to move to a different state (with each new state being equally likely).

This problem is actually fairly straightforward: the stationary distribution for the process x (or y) follows is simply that x is equally likely to be in any of the given states, so P(x = k) = 1/a. Because the initial value for x is also a uniform random distribution, if follows that x is already in the stationary distribution from the start.

It may sound strange, but the value of m has no influence on this. For example, try this code:

StationaryDistribution[
   DiscreteMarkovProcess[
      {1,0,0,0},
      { 
        {1 - m, m/3, m/3, m/3},
        {m/3, 1 - m, m/3, m/3},
        {m/3, m/3, 1 - m, m/3},
        {m/3, m/3, m/3, 1 - m}
      }
   ]
]

The matrix is the transition matrix of the Markov process (for 4 states in this case). Each entry M_ij is the probability that the next state is j if the current state is i (note that the rows of M sum to 1). As you can see by executing the code, the stationary distribution does not depend on m, but is simply 1/4 for every state. This knowledge should make your problem fairly easy to do, I think.

Of course, if the transition matrix is different (for example, not every state is equally likely to be visited), then the story becomes different. I suggest you read up on Markov chains and processes if you want to know more.