{'fraction_non_alphanumeric': 0.04837539305006853, 'fraction_numerical': 0.009003198150984977, 'mean_word_length': 4.105653128430296, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Capacity of DNA Data Embedding Under Substitution Mutations

A number of methods have been proposed over the last decade for encoding information using deoxyribonucleic acid (DNA), giving rise to the emerging area of DNA data embedding. Since a DNA sequence is conceptually equivalent to a sequence of quaternary symbols (bases), DNA data embedding (diversely called DNA watermarking or DNA steganography) can be seen as a digital communications problem where channel errors are tantamount to mutations of DNA bases. Depending on the use of coding or noncoding DNA hosts, which, respectively, denote DNA segments that can or cannot be translated into proteins, DNA data embedding is essentially a problem of communications with or without side information at the encoder. In this paper the Shannon capacity of DNA data embedding is obtained for the case in which DNA sequences are subject to substitution mutations modelled using the Kimura model from molecular evolution studies. Inferences are also drawn with respect to the biological implications of some of the results presented.


I. INTRODUCTION
The last ten years have witnessed the proposal of numerous practical methods [1], [2], [3], [4], [5], [6], [7], [8], [9] for encoding nongenetic information using DNA molecules as a medium both in vitro and in vivo. A conspicuous use of these techniques recently took place when Craig Venter's group produced the first artificial bacteria including "watermarked" information [10]. All of these information encoding proposals hinge on the fact that DNA molecules -which encode genetic information in all living organisms, except for some viruses-are conceptually equivalent to sequences of quaternary symbols.
Therefore DNA data embedding is in essence an instance of digital communications in which channel   ATG TTC CCA AGC ACA TGG TAC GTA TAA   GCC AGG AAT GAT TGT CAG GAG GGC CAT ATC CTC AAG  TTT CCC AGT ACC  TAT GTC TAG   GCT CGA  GGT  ATT CTT  CCT TCA ACT  GTT TGA   GCG CGC  GGG  CTG  CCG TCC ACG  GTG   CGT  TTA  TCT   CGG  TTG   errors are tantamount to mutations of DNA components. The two broad fields of application of DNA data embedding techniques are: 1) the use of DNA strands as self-replicating nano-memories able to store huge amounts of data in an ultra-compact way; and 2) security and tracking applications for genetic material afforded by embedding nongenetic information in DNA (DNA watermarking, steganography, and fingerprinting).
The most basic information theoretical issue in DNA data embedding is the establishment of the upper limit on the amount of information that can be reliably embedded within DNA under a given level of mutations, that is, its Shannon capacity [11]. In this paper we obtain the capacity of DNA data embedding under substitution mutations -which randomly switch the value of bases in a DNA sequencemodelled through a symmetric memoryless channel which was firstly used to study molecular evolution by Kimura [12]. The capacity problem can be straightforwardly tackled when no side information is used by the encoder. The side-informed scenario requires more attention for reasons that will become clear later, and thus occupies us for the best part of this paper. Some biological implications at large of these information theoretical results are also discussed. In particular, the non side-informed scenario happens to be closely connected to previous studies by May, Battail, and other authors that have tried to apply information theoretical concepts to molecular biology.

II. PRELIMINARY CONCEPTS AND ASSUMPTIONS
Chemically, DNA is formed by two backbone strands helicoidally twisted around each other, and mutually attached by means of two base sequences. The four possible bases are the molecules adenine, cytosine, thymine, and guanine, abbreviated A, C, T and G, respectively. Only the pairings A-T and C-G can exist between the two strands, which is why each of the two base sequences is completely determined by the other, and also why the length of a DNA molecule is measured in base pairs (bp). According to this brief description, the interpretation of DNA as a one-dimensional discrete digital signal is straightforward: any of the two strands constitutes a digital sequence formed by symbols from a quaternary alphabet.
As regards the biological meaning of DNA, for the purposes of our analysis it suffices to know that codons -the minimal biological "codewords"-are formed by triplets of consecutive bases in a base sequence. Given any three consecutive bases there is no ambiguity in the codon they stand for, since there is only one direction in which a base sequence can be read. In molecular biology this is called the 5'-3'  Table I, which will be explained in more detail later. The genetic code effectively implements built-in redundancy in terms of protecting protein expression.
The genome of an organism is the ensemble of all its DNA. Segments of a genome that can be translated into proteins through the process described above are called coding DNA (cDNA), whereas those segments that never get translated are called noncoding DNA (ncDNA). A gene is a cDNA segment, or group of segments, which encodes one single protein, and which is flanked by certain start and stop codons (see Table I) plus other markers.
Finally, for each base sequence there are three different reading frames which determine three different codon sequences. The correct reading frame is marked by the position of a start codon.
The main assumptions that we will make in our analysis are the following ones: • ncDNA can be freely appended or overwritten. Although ncDNA does not encode genes, this assumption does not always hold true. This is because certain ncDNA regions act as promoters for gene expression, or are transcribed into regulatory RNA (but not translated into proteins).
However this working hypothesis is valid in suitably chosen ncDNA regions, as proved by several researchers [4], [7] employing live organisms.
• cDNA can be freely modified as long as the genetic code is observed. This is the classic standard assumption supporting the validity of the genetic code. In practice living organisms feature preferred codon statistics, which, if modified, might alter gene expression (for instance translation times, among other effects). Therefore we will also discuss codon statistics preservation in our analysis.
Finally we must mention that we will only consider nonoverlapping genes (either on the same or on opposite strands). However overlapping genes are in any case rare occurrences, except in very compact genomes.
Notation. Calligraphic letters (X ) denote sets; |X | is the cardinality of X . Boldface letters (x) denote row vectors, and 1 is an all-ones vector. If a Roman letter is used both in uppercase (X) and lowercase (x), the two forms denote a random variable and a realisation of it, respectively. p(X = x) is the probability mass function (pmf) of X; we will simply write p(x) when the variable is clear from the context. E[X] is the mathematical expectation of X, and H(X) its entropy. Also, h(q) is the entropy of a Bernoulli(q) random variable. I(X; Y ) is the mutual information between X and Y . Logarithms are base 2, unless explicitly indicated otherwise. The Hamming distance between vectors x and y is denoted by d H (x, y).
A ncDNA sequence will be denoted by a vector x b = [x 1 , x 2 , · · · , x n ], whose elements are consecutive bases from a base sequence. That is, x i ∈ X {A, C, T, G}, the 4-ary set of possible bases. A cDNA sequence will be denoted by a vector of vectors x c = [x 1 , x 2 , · · · , x n ] whose elements are consecutive codons from one of the two antiparallel base sequences, assuming a suitable reading frame among the three possible ones. Therefore, x i ∈ X 3 . We denote by x ′ i α(x i ) ∈ X ′ the amino acid into which a codon x i uniquely translates, which is further discussed below. Also x ′ = α(x c ) = [x ′ 1 , x ′ 2 , · · · , x ′ n ] denotes the unique amino acid sequence established by x c , usually called the primary structure. Using the standard three-letter abbreviations of the amino acid names, we define the set X ′ {Ala, Arg, Asn, Asp, Cys, Gln, Glu, Gly, His, Ile, Leu, Lys, Met, Phe, Pro, Ser, Thr, Trp, Tyr, Val, Stp}. The subset of codons associated with amino acid x ′ ∈ X ′ , that is, S x ′ {x ∈ X 3 |α(x) = x ′ }, is established by the genetic code shown in Table I. The ensemble of stop codons, that is, the stop symbol Stp, is loosely classed as an "amino acid" for notational convenience, although it does not actually map to any compound but rather indicates the end of a gene. We call the number of codons |S x ′ | mapping to amino acid x ′ the multiplicity of x ′ . Due to the uniqueness of the mapping from codons to amino acids, see that S x ′ ∩ S y ′ = ∅ for

A. Mutation channel model
As mentioned in the introduction, an information-carrying DNA molecule undergoing mutations can be readily seen as a digital signal undergoing a noisy communications channel, which we may term "mutation channel" in this context. We will only consider herein substitution mutations (also called point mutations), that is, those that randomly switch letters from the DNA alphabet. We will assume that mutations are mutually independent, which is a worst-case scenario in terms of capacity. Therefore we are assuming that the channel is memoryless and thus accepts a single-letter characterisation; consequently, we will drop vector element subindices whenever this is unambiguous and notationally convenient.
We will model the channel by means of the two-parameter Kimura model of nucleotide substitution [12]. This consists of a 4 × 4 transition probability matrix Π = [p(Z = z|Y = y)], where z, y ∈ X , and which presents the following structure: From this definition, the probability of base substitution mutation, or base substitution mutation rate, is for any y ∈ X , whereas it must hold that 0 ≤ γ ≤ 3/2 so that row probabilities add up to one. The particular structure of Π aims at reflecting the fact that DNA bases belong to one of two categories according their chemical structure: purines, R {A, G}, or pyrimidines, Y {C, T}. There are two types of base substitutions associated to these categories, which in biological nomenclature are: • Base transitions: those that preserve the category which the base belongs to. In this case the model establishes that p(Z = z|Y = y) = (1 − 2γ/3)q for z = y when either both z, y ∈ R or both z, y ∈ Y.
• Base transversions: those that switch the base category. In this case the model establishes that p(Z = z|Y = y) = (γ/3)q for z = y when z ∈ Y and z ∈ R, or vice versa.
The channel model (1) can incorporate any given transition/transversion ratio ε by setting γ = 3/(2(ε + 1)). Estimates of ε given in [13] for the DNA of different organisms range between 0.89 and 18.67, corresponding to γ between 0.07 and 0.79. This range of ε reflects the fact that base transitions are generally much more likely than base transversions due to the chemical similarity among compounds in the same category, that is, ε > 1/2 virtually always in every organism, and therefore γ < 1. However many mutation estimation studies focus only on the determination of q (see for instance [14]), and then one may assume the simplification γ = 1 in the absence of further details. We will make observations at several points for this particular case, which is known as the Jukes-Cantor model in molecular evolution studies. In this situation all off-diagonal entries of Π are equal, that is, p(Z = z|Y = y) = q/3 for all z = y.
Note that the mutation model that we have chosen implies a symmetric channel, since all rows (columns) of Π contain the same four probabilities. Among the memoryless models used in molecular evolution, the Kimura model is the one with higher number of parameters which still yields a symmetric channel. As it is well known, this is advantageous in capacity computations and will be exploited whenever possible.
In the most general case a time-reversible substitution mutations model may have up to 9 independent parameters, and yield a nonsymmetric channel. However according to Li [15] mutation models with many parameters are not necessarily accurate, due to the estimation issues involved.
Under m cascaded mutation stages we have a Markov chain Y → Z (1) → Z (2) → · · · → Z (m) , and model (1) leads to the overall transition probability matrix Π m between Y and Z (m) . As Π = Π T we can and V a matrix whose columns are the normalised eigenvectors of Π associated to the corresponding eigenvalues in D, that is From the diagonalisation of Π m it is straightforward to see that the elements of its diagonal all take the value 1 4 (1 + 2µ m + λ m ), the elements of its skew diagonal take the value 1 4 (1 − 2µ m + λ m ), and the rest of its entries are 1 4 (1 − λ m ). Therefore any row (column) of this matrix contains the same probabilities, as Π m is also the transition matrix of a symmetric channel. From the diagonal elements one can see that the accumulated base substitution mutation rate after m cascaded stages is given by When q > 0, lim m→∞ q (m) | γ>0 = 3/4 but lim m→∞ q (m) | γ=0 = 1/2, because |µ| < 1 for any γ and |λ| < 1 when γ > 0, but λ = 1 when γ = 0. The behaviour of this particular case is connected to the fact that we must have both q ∈ (0, 1] and γ ∈ (0, 3/2] for the Markov chain to be aperiodic and irreducible, and thus possess a limiting stationary distribution. From the previous considerations, the limiting distribution -that is, the distribution of Z (∞) -is uniform, because lim m→∞ Π m = 1 4 1 T 1. When γ = 1 and q = 3/4 then Π = 1 4 1 T 1, and hence every Z (m) is uniformly distributed as in the limiting case.
Lastly, under the base substitution mutation model that we are considering, codons undergo a mutation channel modelled by the 64 × 64 transition probability matrix where ⊗ is the Kronecker product. This is because p(Z = z|Y = y) = 3 i=1 p(Z = z i |Y = y i ) according to our memoryless channel assumption. Trivially this channel is also symmetric. When m mutation stages

A. Noncoding DNA
We will firstly consider this simple case, which will also establish a basic upper bound to cDNA capacity. As per our discussion in Section II, we are assuming that a embedder can overwrite or append a host ncDNA strand x b , which amounts to freely choosing the input y b to the mutation channel. Therefore in this case the channel capacity is given by C nc max I(Z (m) ; Y ) bits/base, where the maximisation is over all distributions of Y . For the mutation model considered, this capacity is that of the symmetric channel, in which H(Z (m) |Y ) is independent of the input and uniformly distributed Y leads to uniformly distributed Z (m) . Hence where for any y ∈ X , that is, the entropy of any row of Π m . Therefore As long as the Markov chain is aperiodic and irreducible then lim m→∞ C nc = 0. The reason is that since the limiting distribution is independent of Y , then lim m→∞ H(Z (m) |Y ) = H(Z (∞) ) = log |X |. It is interesting to note that, under aperiodicity and irreducibility of the Markov chain, this zero limiting capacity will also apply to models more involved than (1), such as those in which the channel matrix is parametrised by up to 9 independent values. Lastly, we also have that C nc | γ=1,q=3/4 = 0, since in this case Z (m) is always uniformly distributed.
As a function of γ the ncDNA capacity is bounded as follows Although it can be shown with some effort that these inequalities always hold true, it is much simpler to prove them for the range of interest γ ≤ 1 and q ≤ 1/2. The latter condition implies that both 0 ≤ λ ≤ 1 and 0 ≤ µ ≤ 1. For fixed m and q, the maximum (respectively, minimum) of C nc over γ corresponds to the minimum (respectively, maximum) of the accumulated base mutation rate q (m) . Differentiating (6) we (as this corresponds to λ ≥ µ), and then its maximum in that range occurs when γ = 1 and its minimum when γ = 0.
The upper bound can be written as Notice that lim m→∞ C nc | γ=0 = 1, that is, the capacity limit is not zero when γ = 0 because then the Markov chain is reducible. This case cannot happen in practice since it would imply that transversion mutations are impossible, but it illustrates that the higher the transition/transversion ratio ε, the higher the capacity.  [14]. We observe the validity of the bounds (10) and the limiting behaviours discussed. From these figures we can also empirically see that a rule-of-thumb capacity cut-off point is given by m ∼ 6/(5γq).

a) Biological interpretations:
We would like to point out that expression (8) also gives the maximum mutual information between a DNA strand and its mutated version in natural scenarios, independent of DNA data embedding procedures. This has sometimes been termed the capacity of the genetic channel in studies applying information theory to molecular biology. Several authors have used the Jukes-Cantor Cascaded mutation stages (m) C nc , bits/base Cascaded mutation stages (m) C nc , bits/base In any case, none of the aforementioned approaches reflects the capacity increase afforded by a mutation model allowing γ < 1. It is possible that the trend towards higher capacity observed as γ → 0 implies that evolution has favoured genetic building blocks which feature an asymmetric behaviour under mutations (in our case, pyrimidines versus purines instead of a hypothetically perfectly symmetric set of four bases for which γ = 1). If this assumption is correct, this symmetry breaking must have occurred early in evolutionary terms, since it is widely believed that the current genetic machinery evolved from a former "RNA world" [21] in which life would only have been based on the self-replicating and catalysing properties of RNA. In the RNA world there would not have been translation to proteins, and therefore no genetic code, and hence information was freely encoded using a 4-ary alphabet almost exactly like the one used in DNA. Note that uracil, which replaces thymine in RNA, is also a pyrimidine, that is, in the RNA world Y = {U, C}. With these facts in mind, we may model the maximum transmissible information under mutations in the RNA world by relying on (8), and thus see that the symmetry breaking conjecture above applies to the evolution of RNA from predecessor genetic building blocks. We must bear in mind that single-stranded molecules, such as RNA, are much more mutation-prone than double-stranded ones such as DNA 1 . Therefore smaller values of m would have sufficed for some type of symmetry breaking to be relevant in terms of information transmission at early stages of life.

B. Coding DNA
Unlike in the ncDNA case, embedding information in cDNA is a problem of coding with side information at the encoder. Given a host sequence x c , the encoder has to modify this host to produce an information-carrying sequence y c which must also encode the same primary structure as x c according to the genetic code. This is equivalent to hiding data in a discrete host under an embedding constraint.
Nevertheless, apart from the trivial difference of using a 4-ary instead of typically a 2-ary alphabet, several issues set apart cDNA data embedding as a special problem. In order to illustrate these issues consider momentarily a typical data hiding scenario in which a discrete binary host, that is x = [x 1 , · · · , x n ] with x i ∈ X = {0, 1}, is modified to embed a message m from a certain alphabet. The watermarked signal y = e(x, m) must be close to x, where closeness is usually measured by means of the Hamming distance d H (y, x). Pradhan et al. [23] and Barron et al. [24] have determined the achievable rate in this scenario, assuming that the elements of X are uniformly distributed, using the average distortion constraint 1 n E[d H (Y, X)] ≤ d, and supposing that y undergoes a memoryless binary symmetric channel with crossover probability q. Their result is where u.c.e{·} is the upper concave envelope. Similarly, our initial goal for cDNA data embedding is obtaining the achievable rate for a fixed distribution of X ′ = α(X) under the symmetric channel discussed in Section II-A, in particular when X is uniformly distributed as in the analyses of Pradhan et al. [23] and Barron et al. [24]. Furthermore we will also obtain capacity, that is, the maximum achievable rate over all distributions of the host X ′ .
The first important difference in the cDNA data embedding scenario is that average inequality constraints on the Hamming distance -such as the ones used in [23], [24]-are meaningless if one wants to carry through to y c the full biological functionality of x c . Instead, since it must always hold that α(y c ) = α(x c ), one must establish the deterministic constraint This requires that d H (y ′ i , x ′ i ) = 0 for all i = 1, · · · , n. The second distinguishing feature of cDNA data embedding is due to the variable support of the channel input variable. Whereas in discrete data hiding with binary host one always has that y i ∈ {0, 1} independently of x i , in cDNA data embedding we have that y i ∈ S α(xi) so that the constraint (12) can always be satisfied. Therefore the support of y i is dependent on x i , as codon equivalence is not evenly spread over the ensemble of amino acids (see Table I).

1) Achievable Rate:
Since side information at the encoder must be taken into account in the cDNA case, then the achievable rate is given by Gel'fand and Pinsker's formula [25] R X ′ c = max I(Z (m) ; U) − I(X ′ ; U) bits/codon, where the maximisation is for nonnegative values of the functional on all distributions p(y, u|x ′ ) under the constraint d H (α(y), x ′ ) = 0, with U an auxiliary random variable that we will discuss next. Note that R X ′ c represents the maximum achievable rate when the host cDNA amino acid sequence is distributed as X ′ .
Gel'fand and Pinsker showed in [25] that in the maximisation problem above one may assume that the channel input is a deterministic function of the side information X ′ and the auxiliary variable U, that is, Y = e(X ′ , U). Since the support of Y|x ′ must be the set of codons S x ′ corresponding to amino acid x ′ -so that the biological constraint can always be satisfied-then the cardinality of the support of U|x ′ has to coincide with the multiplicity of x ′ , that is, |S x ′ |. The support of U|x ′ must actually be S x ′ , because U must also act as a good source code for X ′ in order to minimise I(X ′ ; U) under the genetic constraint, and if the support of U|x ′ is otherwise then the constraint cannot always be met. One can now establish Y|x ′ = U|x ′ without loss of generality, although any permutation of the elements of S x ′ is actually valid to define Y|x ′ = e(x ′ , U). Therefore in the following one may consider that Y = U, that is, that U is the mutation channel input. Noticing that S x ′ ∩ S y ′ = ∅ for x ′ = y ′ ∈ X ′ , the distribution of U can be put as p(u) = p(u|x ′ )p(x ′ ) when u ∈ S x ′ . This discussion on U also implies that H(X ′ |U) = 0, since given a codon u there is no uncertainty on the amino acid represented, and therefore I(X ′ ; U) = H(X ′ ) Since Y|(x ′ , u) is deterministic, from the considerations above we have that the achievable rate for a fixed distribution of X ′ is given by As H(Z (m) |U) only depends on the transition probabilities of the symmetric channel, and as trivially H(X ′ ) only depends on X ′ , (13) amounts to the constrained maximisation of H(Z (m) ).
There are several cases in which (13) can be analytically determined, which are discussed next. First of all, since C nc | γ=1,q=3/4 = 0 then R X ′ c | γ=1,q=3/4 = 0 for any X ′ , because R X ′ c ≤ 3C nc . Therefore in this catastrophic case the choice of p(u|x ′ ) is irrelevant. Furthermore it can be shown that p(u|x ′ ) = 1/|S x ′ |, that is, U|x ′ uniformly distributed, is the maximising strategy in two situations, which are discussed in the following lemmas.

Lemma 1. If q = 0 then the achievable rate is
Proof: Using the chain rule of the entropy we can write H(U, As H(X ′ |U) = 0, and as Z (m) = U when q = 0, then the achievable rate is given by We just need to see now that H(U|X ′ ) = is maximum for all x ′ , which implies that U|x ′ be uniformly distributed in all cases. Then H(U|x ′ ) = log |S x ′ | and (14) follows. (14) is the embedding rate intuitively expected in the mutation-free case. For example, if X were uniformly distributed, which would yield X ′ = α(X) nonuniform with pmf p(x ′ ) = |S x ′ |/|X | 3 , then we would obviously compute the rate as R

Remark. Note that
7819 bits/codon, since |S x ′ | choices are available to the embedder when the host amino acid is x ′ . The rate in the uniform case can actually be obtained in closed form for every q using the following result.

Lemma 2. If X is uniformly distributed then the achievable rate is
where C nc max I(Z (m) ; U) and this maximisation is unconstrained on p(u), that is, C nc is the capacity of the symmetric codon mutation channel.
Proof: Since p(u) = p(u|x ′ )p(x ′ ) when u ∈ S x ′ , with uniformly distributed X we have that p(u) = p(u|x ′ )|S x ′ |/|X | 3 when u ∈ S x ′ . Therefore choosing U|x ′ to be uniformly distributed implies that p(u) = 1/|X | 3 for all u. Since Π m is symmetric and a uniform input maximises mutual information over a symmetric channel, then C nc = max I(Z (m) ; U) is achieved in (13).
Remarks. Since C nc | q=0 = log |X | 3 , observe that the particular case in the previous remark can be written as well as R An interesting insight is also afforded by seeing that the three parallel symmetric channels undergone by the bases in a codon are mutually independent, and hence one can use the equality C nc = 3 C nc in (15). As H(X ′ ) is the lower bound to the lossless source coding rate of X ′ , expression (15) tells a fact that is intuitively appealing but which is only exact when X is uniform: the cDNA embedding rate is the same as three times the ncDNA embedding rate minus the rate needed to losslessly convey the primary structure of the host to the decoder.
Unlike in the case considered above, the distribution of X in real cDNA sequences (that is, genes) is not uniform. To start with, there can only be a single Stp codon in a sequence that encodes a protein (gene). As in many other channel capacity problems, it does not seem possible in general to analytically derive the optimum set of pmf's p(u|x ′ ) in order to compute the achievable rate R X ′ c corresponding to a host distributed as X ′ . To see why one can pose the analytical optimisation problem and see that it involves solving a nontrivial system of |X | 3 + |X ′ | nonlinear equations and unknowns. However the numerical solution is straightforward by means of the Blahut-Arimoto algorithm [26] adapted to the side-informed scenario. Such an algorithm has been described by Dupuis et al. [27]. An example of the optimal distributions numerically obtained for a particular case is shown in Figure 3.
A last observation is that, in general, p(u|x ′ ) = 1/|S x ′ | turns out to yield a good approximation to the exact numerical solution. Note that for distributions of X ′ different from the one in Lemma 2 one cannot produce a uniform input U to the symmetric channel in order to generate a uniform output Z (m) . This is the case illustrated in Figure 3; nonetheless, note from this figure that p(u|x ′ ) does not differ excessively from a uniform distribution for several amino acids x ′ . A justification of this behaviour is as follows. Using the fact that conditioning cannot increase entropy, a suboptimal maximisation approach is given by maximising the lower bound ). This requires maximising H(Z (m) |x ′ ) = − z∈X 3 p(z|x ′ ) log p(z|x ′ ) for all x ′ . Observing from Table I that codons mapping to the same amino acid share in many cases up to two bases, we can approximate With this approximation, whenever u∈S x ′ p(z|u) is constant for all z ∈ S x ′ , choosing p(u|x ′ ) to be uniform implies that p(z|x ′ ) is also uniform, which maximises H(Z (m) |x ′ ). It can be verified that this condition holds for all x ′ such that |S x ′ | = 1, 2, 4, which accounts for 16 out of the 21 elements in X ′ .  GCG  AGA  AGG  CGA  CGC  CGT  CGG  AAC  AAT  GAC  GAT  TGC  TGT  CAA  CAG  GAA  GAG  GGA  GGC  GGT  GGG  CAC  CAT  ATA  ATC  ATT  CTA  CTC  CTT  CTG  TTA  TTG  AAA  AAG  ATG  TTC  TTT  CCA  CCC  CCT  CCG  AGC  AGT  TCA  TCC  TCT  TCG  ACA  ACC  ACT  ACG  TGG  TAC  TAT  GTA  GTC  GTT  GTG  TAA  TAG   Conditional pmf's are depicted in alternating red and blue colours to facilitate plot reading. distributions corresponding to two real genes: Ypt7 (S. Cerevisiae) and FtsZ (B. Subtilis), whose GenBank accession numbers are NC 001145 and NC 000964, respectively; also depicted are the rate (15) for X uniform and the rate for the deterministic distribution of X ′ with outcome Ser, which, as we will discuss in Section III-B2, yields capacity. We observe in these plots that there is barely a difference in the results obtained with the Blahut-Arimoto algorithm and the uniform approximation to p(u|x ′ ) that we have discussed.
a) Codon statistics preservation: If we require that the original codon statistics of the host are preserved in the information-carrying sequence, then we must peg p(u|x ′ ) to the corresponding distribution of the host. Therefore in this case no maximisation is required, and the corresponding rate will be lower or equal than the one achieved without codon statistics preservation. This type of constraint is equivalent to a steganographic constraint in data hiding, since the pmf of the host is preserved in the information-carrying sequence. A comparison of maximum rates and codon statistics preservation rates for the same genes as before is given in Figure 8. Note that in the uniform case both rates coincide because of Lemma 2. 2) Capacity: It remains the computation of capacity, that is, It is simple to explicitly obtain C c in two particular cases discussed in Section III-B1.
• γ = 1 and q = 3/4: obviously, C c | γ=1,q=3/4 = R X ′ c | γ=1,q=3/4 = 0. However, the point that has to be made here is that, although this is true for any X ′ , only deterministic X ′ yields H(X ′ ) = 0 exactly, and hence, by the continuity of the rate functional, this will be the best strategy when approaching q = 3/4 from the left.

Remark. A trivial upper bound for any
side-informed cDNA data embedding capacity will not be able to achieve non-side-informed ncDNA capacity for every mutation rate. This is similar to parallel results in side-informed encoding with discrete hosts [23], [24] (for uniform side information), and unlike the well-known result by Costa for continuous Gaussian hosts [28].
From our previous discussion on the value of C c for two particular cases one may conjecture that a pmf with support in W ′ may be capacity-achieving. The actual capacity-achieving strategy is given by the following theorem: Theorem 1. Capacity is achieved by the deterministic pmf of X ′ that maximises H(Z (m) ).

a) Biological interpretations:
The results in this section solely concern artificial embedding of information in cDNA, and thus seem to have less obvious applicability in biological terms than the ones concerning ncDNA. However an intriguing phenomenon which somehow indicates a biological connection of these results can be observed in Figures 11-12 which concern achievable rates for sequences encoding a single amino acid/symbol. The effect is observed for γ < 1 -that is, the range of γ in which the model is more realistic-and consists of a rate droop for two particular values of ξ ′ as m → ∞ with respect to all other symbols presenting the same multiplicity. The particularity is that these two values of ξ ′ correspond to the stop symbol (Stp) and the amino acid Leu, which happens to double as start codon in prokaryotes. Therefore all stop codons and most of the start codons seem to be less suited to carrying extra information (redundancy) when isolated. It is not obvious how to interpret this effect, but one may surmise that these codons might have suffered some type of selective pressure during the emergence of the genetic code which somehow depended on the information theoretical amount studied here. In the case of the stop symbol, this may be due to the fact that it can only appear once per gene, and so it makes for a bad conveyor of extra information beyond its basic function.

IV. CONCLUSIONS
We have provided an analysis of the embedding capacity of DNA when mutations are modelled according to the Kimura model from molecular evolution studies, and discussed some biological connections of these results. A more thorough study would require considering insertion and deletion mutations (indels). Although the exact computation of capacity under indels is an unsolved problem in most digital communications scenarios, some approximations relying on realignment methods from bioinformatics might suffice in this context. Generalisations of the Kimura model may also be considered. Although in general they will lead to nonsymmetric channels, these can be numerically handled using the Blahut-Arimoto algorithm.

APPENDIX
A. Capacity-achieving strategy p(x ′ ) In order to find the capacity-achieving strategy we need to solve for x ′ ∈ X ′ , with ν a Lagrange multiplier. In the following we will write p(z|x ′ ) = p(Z (m) = z|X ′ = x ′ ) for notational convenience. Assuming natural logarithms for simplicity, and using ∂p(z)/∂p(x ′ ) = p(z|x ′ ), (17) becomes for x ′ ∈ X ′ . The solution remains unchanged if we multiply (18) across by p(x ′ ). This allows us to see by inspection that any extreme of the Lagrangian in (17) has to be deterministic, that is, p(x ′ ) = 1 for some x ′ = ξ ′ and p(x ′ ) = 0 for x ′ = ξ ′ . Note that this is in agreement with the strategies for the cases q = 0 and γ = 1 with q = 3/4 discussed in Section III-B2. See for instance that a uniform distribution of X ′ cannot possibly solve (18) for all x ′ ∈ X ′ , because z p(z|x ′ ) log y ′ p(z|y ′ ) is not constant on x ′ unless γ = 1 and q = 3/4, in which case we have shown that capacity is zero for any distribution.
According to the previous discussion, for any capacity-achieving solution it always holds that H(X ′ ) = 0, and then we just have to maximise H(Z (m) ) over the ensemble of 21 deterministic distributions of The computation of R ξ ′ c and of the maximising distribution U|ξ ′ can be done using the Blahut-Arimoto algorithm, following the discussion in Section III-B1 on the optimal strategy for fixed p(x ′ ). Note that ξ ′ = Trp and ξ ′ = Met can be ruled out outright, since |S Trp | = |S Met | = 1, and then only null rates are possible in these cases. Then we only need to compute R ξ ′ c for 19 amino acids. Also, ξ ′ = Stp can only be considered hypothetically, since this symbol can only appear exactly once in a gene.
a) Approximation to maximising strategy: It is also possible to provide a closed-form approximation to the maximising distribution U|ξ ′ , which yields a better approximation to the embedding rate than just using the approximation p(u|ξ ′ ) = 1/|S ξ ′ | discussed in Section III-B1. Observe firstly that when X ′ is deterministic the situation is equivalent to a non-side informed discrete channel with |S ξ ′ | inputs and |X | 3 outputs, with a transition probability matrix Λ whose rows are the rows of Π m corresponding to the codons associated with ξ ′ . In general this channel will not be symmetric nor weakly symmetric, since although its rows are permutations of the same set of probabilities, its columns are not, and their sum is not constant either. However H(Z (m) |U) is still independent of the distribution of U, and then we only need to maximise H(Z (m) ) to find capacity. The corresponding conditions for the maximum are z∈X 3 p(z|v) log p(z) + 1 = ρ, for v ∈ S ξ ′ , and with ρ a Lagrange multiplier.
This yields the linear system where the probabilities p(u|ξ ′ ), with u ∈ S ξ ′ , are the elements of the 1 × |S ξ ′ | vector π (arranged in the same codon order as the rows of Λ), and 1 is an all-ones vector of size 1 × |S ξ ′ |. Since π must be a pmf, we may fix any arbitrary value of ρ, such as ρ = 1, and then normalise the solution π to the resulting linear system, that is The matrix ΛΛ T is invertible if both q = 1/(4γ/3) and q = 1/(2(1 − γ/3)) because in this case the rows of Λ are linearly independent. This is due to the fact that under the two conditions above the rows of Π m are linearly independent, since its eigenvalues are all the possible products of three eigenvalues of Π m [16] and the conditions above guarantee that these are nonzero. A sufficient condition for the invertibility of ΛΛ T is q < 1/2, which spans most cases of interest.
Since we have linearised the optimisation problem then π may contain negative values, but in practice these are relatively small. Setting these values to zero and normalising π we obtain an approximation to the optimum distribution p(u|ξ ′ ). An example of this approximation compared to the results of the Blahut-Arimoto algorithm is shown in Figure 13.