text:
# Many-body interactions and correlations in coarse-grained descriptions of polymer solutions

## I Introduction

An efficient statistical description of condensed matter systems and materials almost invariably involves some degree of coarse-graining, whereby a large fraction of the initial microscopic degrees of freedom are traced out, leaving a much reduced space of variables associated with the composite entities or pseudo-particles. The latter are then coupled via effective interactions which result from the partial averaging over the initial microscopic degrees of freedom. The reduction of the initial multicomponent system to a coarse-grained system with a substantially smaller set of composite particles implies that the resulting effective interactions may involve three-body and higher-order contributions, even if the original multicomponent system involved only pair-wise additive forces, like Coulombic interactions. Alternatively, one may wish to retain the simplicity of pair-wise additivity of the effective interactions, but the price to pay is that such effective pair potentials are then state-dependent, e.g. are functions of the temperature and/or density. The reason for this is that the effective interaction energy is a free energy associated with the averaged-out degrees of freedom, which generally has an entropic component.

There are many examples of the coarse-graining procedure which has just been outlined. In molecular systems the forces between nuclei result from gradients of the electronic ground state energy surface which depends parametrically on the nuclear coordinates, and adjusts adiabatically to the slow motion of the latter within the Born-Oppenheimer approximation. This scenario is mimicked, at least at the level of valence electrons, in “ab initio” Molecular Dynamics simulations pioneered by Car and Parrinello. However in situations where no strong covalent or hydrogen bonding is present, the more common route is to represent the total ground state electronic energy surface by a sum of one-body, two-body and higher order terms. The one-body contribution is the sum of the ground state energies associated with individual, isolated chemical entities (atoms, ions, or molecules). The two-body term is made up of the sum of pair potentials acting between molecules and higher order terms correspond to isolated clusters of three or more molecules. The sum of ground state energies of individual molecules does not contribute to the forces between them, and can hence be ignored in the description of collective equilibrium or transport processes which do not involve chemical reactions. In the simplest case of rare-gas atoms, the pair-wise interactions would include overlap repulsion at short range and dispersion forces at long range, while triplet interactions would include, among others, the Axilrod-Teller triple dipole dispersion potential, which contributes very significantly to the thermodynamic and transport properties of the heavier rare gases in their condensed states. The effect of the higher order interactions can be approximately incorporated into an effective pair potential, which differs from the bare pair potential, valid for an isolated pair of molecules, and becomes density-dependent.

Similarly in metals an effective interaction between ions may be determined by tracing out the conduction electrons, using perturbation theory or response theory. Treating the ion-electron coupling to lowest order (linear response) leads to a structure-independent volume term and to a pair-wise screened effective potential between “dressed” ions or pseudo-atoms, which both depend on the macroscopic conduction electron density. The two-body level is generally sufficient for alkali and other simple metals, in part because of a quantum interference effect which strongly decreases the magnitude of the higher order response terms. For multivalent and transition metals many-body effective interactions can no longer be neglected, and a full “ab initio” treatment may be necessary.

Coarse-graining becomes crucial in the highly asymmetric systems of soft matter, involving macromolecules or colloidal particles, as well as molecular scale entities, like solvent molecules or ions. The latter are traced out to derive effective interactions between the electric double-layers associated with charged surfaces (colloids, membranes etc…). The microscopic ions play a role similar to valence electrons in metals, but with the important difference that quantum degeneracy effects are absent, and that finite temperature entropic effects control the width of the electric double-layers, with a resulting density and temperature dependence of the effective interactions between the mesoscopic colloidal particles. A good example of a such a pair potential is provided by the classic Derjaguin-Landau-Verwey-Overbeek (DLVO) effective pair potential between spherical charged colloidal particles. Three-body interactions can be derived in similar fashion. However, the averaging over the microion degrees of freedom also leads to a structure-independent, but state-dependent volume term, part of which is associated with the self-energy of individual double-layers, and which is reminiscent of the volume term in metals; this term has a profound effect on the phase diagram of charge-stabilized colloids.

Another important class of effective interactions of entropic origin, which follow from averaging over the configurations of non-adsorbing polymers and small colloidal particles, are the depletion forces which have received much renewed interest in recent years. Depletion pair potentials depend strongly on the concentration of the depletant; recent attempts have been made to compute the three-body interactions from simulations or density functional theory. Volume terms arising from the depletion potential picture may have an important effect on the osmotic equation of state, but they are not expected to influence the phase behavior in these uncharged systems.

The present paper is concerned with a coarse-grained description of dilute and semi-dilute solutions of polymers in good solvent. Whereas tracing out the microscopic ions in a charged colloidal suspension has many analogies with the liquid metal problem, the coarse-graining of neutral polymers, achieved by integrating out the internal monomeric degrees of freedom, resembles more closely the case of effective potentials between neutral atoms and molecules, obtained by tracing out the internal electronic degrees of freedom. The basic idea explored here, which goes back to Flory and Krigbaum, is to represent a set of polymers, each made up of $`L`$ monomers or segments, as single particles, interacting with each other through an effective interaction between their centers of mass (CM). The important point, realized by Grosberg et al. , is that the effective pair interaction remains finite, even for infinitely long polymers. Monte Carlo simulations and renormalization group calculations show that for two isolated non-intersecting polymer coils, the effective potential between their CM’s is of order $`2k_BT`$ in the scaling limit, i.e. for $`L`$ going to infinity, while the range of the interaction is of the order of the radius of gyration, $`R_g`$, of the polymers. Recently we have extended this investigation by simulating large systems of self-avoiding walk (SAW) polymers at finite concentration. The resulting CM pair distribution function $`g(r)`$ was then inverted to yield a concentration-dependent effective pair potential $`v(r;\rho )`$, where $`\rho `$ is the number of polymer coils per unit volume. Although $`v(r;\rho )`$ was not found to change dramatically with $`\rho `$, the $`\rho `$-dependence of $`v`$ is very significant for the accurate determination of the osmotic properties of dilute and semi-dilute polymer solutions.

In this paper we adopt a somewhat different point of view, by determining state-independent effective pair, triplet, quadruplet and quintuplet interactions; these many-body interactions are determined by successively considering clusters of 2, 3, 4, and 5 SAW polymer chains, determining the corresponding n-body distribution functions, from which an effective n-body potential is derived. Triplet interactions between the cores of star-polymers have recently been determined in a similar way. The next step taken in this paper is to relate the low-density (state-independent) pair and triplet interactions to the density-dependent effective interactions determined in our earlier work. A final section will be devoted to an analysis of three-body correlations as measured by bond-angle distributions and deviations from the Kirkwood superposition approximation. The convergence of the series of n-body interactions is assessed on the basis of scaling arguments in the Appendix.

## II Simulation Models and Methods

Many properties of polymers in a good solvent are well described by models which ignore all microscopic details of the intermolecular interactions, except their excluded volume. For that reason polymers are often modeled as self avoiding walks (SAW) on a lattice, a model lending itself well to efficient computer simulations. We consider the situation of $`N`$ athermal SAW chains of length $`L`$ on a simple cubic lattice of $`M`$ sites. The bead or segment concentration is given by $`c=NL/M`$, while the polymer chain concentration is given by $`\rho =N/M`$. The polymers are characterized by their radius of gyration $`R_g`$ which, for an isolated polymer, scales as $`R_gL^\nu `$, where $`\nu 0.59`$ is the Flory exponent. We also define an overlap concentration $`\rho ^{}=1/\frac{4}{3}\pi R_g^3`$ at which there is on average one polymer per sphere of radius $`R_g`$. Solutions with $`\rho /\rho ^{}<1`$ are called dilute, while solutions with $`\rho /\rho ^{}>1`$ and $`c1`$ are called semi-dilute. When the monomer density $`c`$ becomes appreciable, the solution moves from the semi-dilute to the melt regime. In this paper we will focus on densities $`\rho /\rho ^{}2`$, i.e. the dilute regime and the beginning of the semi-dilute regime.

When modelling the semi-dilute regime, it is important to take sufficiently long polymer chains. The first reason is that for studying the semi-dilute regime one needs a large polymer density $`\rho `$ together with a low monomer density $`c`$. We found earlier that the monomer density $`c^{}`$ at the overlap concentration $`\rho ^{}`$ scales roughly like $`c^{}4L^{0.8}`$ for SAW polymers on a simple cubic lattice. Throughout this paper we use polymers of length $`L=500`$ for which $`c0.05`$ at $`\rho /\rho ^{}=2`$, so that we are still clearly in the semi-dilute regime. In contrast, for $`L=100`$ the monomer density is $`c0.2`$ at $`\rho /\rho ^{}=2`$, suggesting that a meaningful semi-dilute regime does not exist for such short polymers.

The second reason for using long polymer chains is that we want to study properties — particularly the effective potentials between polymer chains — in the scaling regime, where all length dependence is completely captured by $`R_g`$. In a previous paper we established that two properties of the effective potential relevant to thermodynamics, namely the second-virial coefficient $`B_2`$ between two polymers and the effective pair-potential between CM’s, multiplied by the square of the CM distance, $`r^2v(r)`$, are very close to the scaling limit for $`L=500`$ polymers, the length we will use in this paper.

Simulations were done with the Monte Carlo pivot algorithm  combined with simple translational moves. For concentrations $`\rho >\rho ^{}`$ we also use Configurational Bias Monte Carlo algorithms . For $`L=500`$ polymers we find that the radius of gyration of an isolated coil is $`R_g=16.50\pm 0.03`$. We used a simulation box of size $`M=(240)^3`$, and varied the number of polymers from $`N=2`$ for the 2-body calculations to $`N=6400`$, which corresponds to $`\rho /\rho ^{}=8.9`$.

Since we are dealing with athermal chains consisting of monomers interacting only via hard-core repulsion, we set the reciprocal temperature $`\beta =1/k_BT=1`$ throughout this paper.

## III Density independent many-body interactions

### A Expanding the coarse-grained free energy in a series of many-body interactions

Following the discussion in, the Helmholtz free energy $``$ of a set of $`N`$ polymers of length $`L`$ with their centers of mass (CM) distributed according to the set of coordinates $`\{𝐫_i\}`$, in a volume $`V`$, can be written as the following expansion:

$`(N,V,\{𝐫_i\})`$ $`=`$ $`^{(0)}(N,V)+{\displaystyle \underset{i_1<i_2}{\overset{N}{}}}w^{(2)}(𝐫_{i_1},𝐫_{i_2},)`$ (1)
$`+`$ $`{\displaystyle \underset{i_1<i_2<i_3}{\overset{N}{}}}w^{(3)}(𝐫_{i_1},𝐫_{i_2},𝐫_{i_3})+\mathrm{}`$ (2)
$`\mathrm{}`$ $`+`$ $`{\displaystyle \underset{i_1<\mathrm{}<i_n}{\overset{N}{}}}w^{(N)}(𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_N})`$ (3)
In the scaling limit, each term in the series is independent of $`L`$ as long as the n-tuple CM coordinates $`\{𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_n}\}`$ are expressed in units of $`R_g`$, the radius of gyration at zero density. Note that this coarse-grained free energy includes an implicit statistical average over all the monomeric degrees of freedom for a fixed configuration $`\{𝐫_i\}`$ of the CM. The full free energy of the underlying polymer system can be calculated as follows:

$$F(N,V)=\mathrm{ln}\underset{\{𝐫_i\}}{}\mathrm{exp}\left[(N,V,\{𝐫_i\})\right]$$
(4)
so that Eq. (1) can be viewed as an expansion of the effective interaction between the CM in terms of (entropic) many-body interactions. $`^{(0)}(N,V)`$ is the so-called volume term, the contribution to the free energy which is independent of the configuration $`\{𝐫_i\}`$. Here it includes the free energy of a single isolated polymer which is independent of the position of its CM in a homogeneous solution; translational invariance also implies that there is no one-body term in the expansion. Each subsequent term $`w^{(n)}(𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_n})`$ is defined as the free energy of $`n`$ polymers with their CM positions at $`\{𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_n}\}`$ minus the contributions of all lower order terms. In other words, it is the contribution to the free energy of $`n`$ polymers which is not included in the sum of all lower order terms. For instance, the 2-body term $`w^{(2)}(r_{ij})`$ can be defined as the difference between the coarse-grained free energy $``$ for two particles with their CM distance held at $`r_{ij}=|𝐫_𝐢𝐫_𝐣|`$, and the free energy of the same two polymers when they are infinitely far apart. Here we use the translational and rotational invariance of a homogeneous system to reduce the number of degrees of freedom. Similarly, the three-body term for a given triplet configuration $`\{𝐫_i,𝐫_j,𝐫_k\}`$ can be written in terms of only 3 variables (see Fig. 1) as

$`w^{(3)}(r_{ij},r_{jk},r_{ki})`$ $`=`$ $`(N=3,V,r_{ij},r_{jk},r_{ki})`$ (5)
$`^{(0)}(N=3,V)`$ $``$ $`w^{(2)}(r_{ij})w^{(2)}(r_{ik})w^{(2)}(r_{jk}).`$ (6)
In other words, it is that part of the effective interaction between three polymers which cannot be described by volume and pair interaction terms alone. In principle, this procedure may be continued until, for a system with $`N`$ polymers, the $`N`$th term determines the total coarse-grained free energy. In practice, this approach is not feasible because the number of $`n`$-tuple coordinates increases rapidly with $`n`$, as does the complexity of each higher order term, so that the series in Eq. (1) quickly becomes intractable. Instead, one hopes to show that the series converges fast enough that only a few low order terms are needed to obtain a desired accuracy. We now turn to the derivation of these density independent potentials for our system of SAW polymers.

### B 2-body interactions

There is a general relationship between the $`\rho 0`$ limit of $`n`$-body correlation functions and the $`n`$-body potential. For the 2-body case this reduces to

$$\underset{\rho 0}{lim}g_2(r)=\mathrm{exp}\left[w^{(2)}(r)\right],$$
(7)
where $`g_2(r)`$ is the pair distribution function. Although this definition resembles that of the potential of mean force (PMF), usually defined as $`w^{PMF}(r)=\mathrm{ln}[g_2(r)]`$ for any density, they are only equivalent in the limit of infinite dilution. Strictly speaking, the PMF is not a potential but simply a restatement of the pair-correlations. Consider the simplest case, namely a system with no higher-order $`(n>2)`$ interactions. If one were to use a finite density PMF as a pair potential at that same finite density, the resulting pair correlations would not be those of the system used to derive the PMF. In contrast, the potential defined in Eq. (7) is the correct pair potential which would exactly reproduce the pair correlations of the original system.

In our simulations we calculate $`w^{(2)}(r)`$ from the logarithm of the overlap probability as a function of CM distance (c.f. Eq. 7). Although the arguments above were made for a free energy in a continuous space, they easily carry over for the lattice model we simulate. In fact, the CM lives on a finer grid than the original SAW lattice polymers, so that our results are already closer to the continuum limit. The overlap probability is determined by sampling the configurations of two polymers infinitely far apart with the pivot algorithm, and, after every 1000 pivot moves, searching for any monomer overlaps as a function of the CM distance. The effective potential calculated in this manner has a near Gaussian shape with a value at full overlap of $`w^{(2)}(0)=1.88\pm 0.01`$ for our $`L=500`$ polymers, very close to the scaling limit estimate $`w^{(2)}(0)=1.80\pm 0.05`$, and a range of the order of $`R_g`$, as shown in Fig. 2. This picture agrees with earlier RG and simulation studies. Note that in the scaling limit the potentials depend only on $`R_g`$, so that the free energy cost of completely overlapping the CM of two polymers is independent of their length $`L`$. That this free energy cost at complete overlap should depend weakly on polymer length follows from their fractal nature, but more sophisticated scaling theory arguments are needed to prove that $`w^{(2)}(0)L^0`$.

### C 3-body interactions

Just as the 2-body interactions follow from the low density limit of the pair correlations, so also the 3-body or triplet interactions between the CM of three polymers can be derived by taking the low-density limit of the 3-body distribution function

$$g_3(𝐫_1,𝐫_2,𝐫_3)=\frac{1}{Z}\mathrm{}e^{(\{𝐫_i\})}𝑑𝐫_4\mathrm{}.d𝐫_N,$$
(8)
where $`(\{𝐫_i\})`$ is the coarse-grained free energy (1) and and $`Z=\mathrm{exp}[F]`$ is the configurational integral for the CM coordinates defined in Eq. (4) the effective interactions already include an average over monomeric degrees of freedom, this configurational integral is equivalent to that of the full polymer system.) Taking the low-density limit for a homogeneous system gives

$$w^{(3)}(r_{12},r_{23},r_{13})=\underset{\rho 0}{lim}\mathrm{ln}\left[\frac{g_3(r_{12},r_{23},r_{13})}{g_2(r_{12})g_2(r_{23})g_2(r_{13})}\right]$$
(9)
The $`g_2(r)`$ ensure that the contributions due to the pair interactions are subtracted from the triplet interaction (c.f. Eq. 5). For our homogeneous system the 3-body potential depends only on the three variables $`\{r_{12},r_{23},r_{13}\}`$ shown in Fig. 1. Even then, calculating the triplet interaction $`w^{(3)}(r_{12},r_{23},r_{13})`$ for every possible triplet arrangement is very cumbersome. We therefore confine ourselves to configurations that make up an equilateral triangle. Instead of three variables, the potential now depends only on the length $`r`$ of each side of the triangle, simplifying the calculation of Eq. (9) to:

$$w^{(3)}(r)=\underset{\rho 0}{lim}\left[\mathrm{ln}g_3(r)3w^{(2)}(r)\right],$$
(10)
where we also used Eq. (7). We expect that for $`r=0`$, i.e. complete overlap, the three body interaction will be strongest, while for large $`r`$ the interaction should vanish.

We calculated $`w^{(3)}(r)`$ for three $`L=500`$ SAW polymers on a lattice. At this infinite dilution, $`g_3(r)`$ is simply the probability that three polymers in a configuration $`\{r_{12},r_{23},r_{31}\}`$ do not overlap. In the Monte Carlo simulation we integrate over the monomeric degrees of freedom by performing pivot moves. Once every 1000 MC steps we move the polymers into a triangular configuration and check for overlap. The results are plotted in Fig. 2. Since the total free energy increases with the number of polymers, a more relevant measure of the 3-body interactions is the relative potential $`\mathrm{\Delta }w^{(3)}(r)=w^{(3)}(r)/3w^{(2)}(r)`$ which denotes the strength of the 3-body interaction relative to that of the two body interactions. As shown in the inset of Fig. 2, the relative contribution of the 3-body potential to the total free energy is quite small, only about $`9\%`$ of the contribution from the pair potentials.

### D 4-body interactions

Next we turn to the 4-body interactions. Again, even for a homogeneous phase, the total number of relevant coordinates makes the calculation of the full interaction prohibitively complex. To restrict the number of coordinates in our calculations we determine the 4-body potential by placing the four polymer CM’s on a regular tetrahedron and determining the non-overlap probability as a function of the length $`r`$ of each side of the tetrahedron. The 4-body potential is then defined as

$$w^{(4)}(r)=\underset{\rho 0}{lim}\mathrm{ln}\left[g_4(r)\right]4w^{(3)}(r)6w^{(2)}(r)$$
(11)
since the tetrahedral configuration includes four equilateral triangles (3-body interactions) and six edges (2-body interactions). The full and relative 4-body interaction, $`\mathrm{\Delta }w^{(4)}(r)=w^{(4)}(r)/(4w^{(3)}(r)+6w^{(2)}(r))`$, are plotted in Fig. 2. Note that the 4-body interaction is smaller in absolute magnitude and has the opposite sign to the 3-body interaction. The relative contribution of the 4-body interactions to the total free energy is less than $`5\%`$ of the total potential, and also less than the relative 3-body contribution.

### E 5-body interactions

Calculating the 5-body interaction is even more complicated than the 4-body interaction, and so we only evaluate it at full overlap – when all the CM’s coincide – where we expect its contribution to be largest. More generally, for any $`n`$th order term the interaction at full overlap is given by:

$$w^{(n)}(0)=\underset{\rho 0}{lim}\left[\mathrm{ln}g_n(0)\right]\underset{m=2}{\overset{m=n1}{}}\left(\genfrac{}{}{0pt}{}{n}{m}\right)w^{(m)}(0)$$
(12)
where $`lim_{\rho 0}g_n(0)`$ is the normalized probability of full overlap of the CM of $`n`$ polymers. As long as the particles are equidistant from each other, the same combinatorial expression holds for finite $`r`$. For the 5-body term we find that $`w^{(5)}(0)=0.4\pm 0.15`$, while the relative contribution of the 5-body terms is given by $`w^{(5)}(0)/(5w^{(4)}(0)+10w^{(3)}(0)+10w^{(2)}(0))=0.027\pm 0.01`$. Again, the relative contribution of the 5-body term to the free energy is smaller and of opposite sign to those of the 4-body terms. Going beyond the 5-body interaction, even at complete overlap becomes increasingly difficult. For example, for the 5-body interaction, of the $`10^8`$ overlap checks, each attempted after $`1000`$ pivot moves, only about $`30`$ resulted in non-overlap. For the 6-body interaction we estimate that $`10^{11}`$ MC attempts would be needed. Another problem arises from finite monomer density. As more and more polymers overlap, the monomer density increases, so that in practice for a given polymer length $`L`$, only a finite number of multiple overlaps are possible. We found previously that the largest finite-size corrections to the scaling limit were at $`r=0`$ for $`w^{(2)}(0)`$. The same probably holds for the higher order interactions. However, our $`L=500`$ calculations should still be very near the scaling limit.

## IV Effective density dependent pair interactions

From the previous section we see that explicitly calculating the density independent interactions $`w^{(n)}(𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_N})`$ becomes rapidly more complex with increasing order $`n`$. Calculating all higher order terms is therefore impossible. In this section, we describe a way to include the average effect of all higher order terms by extending the relationship between the pair-interactions and the pair-correlations to finite density $`\rho `$. This leads to a density dependent effective pair interaction $`v(r;\rho )`$.

### A Inverting pair-correlations to derive density dependent pair potentials

Although for finite densities there is no known direct functional relationship of the type of Eq. (7), there is a theorem which states that for any given pair-correlation function $`g_2(r)`$ and density $`\rho `$, there exists (up to an additive constant) a unique pair potential $`v(r;\rho )`$ which exactly reproduces $`g_2(r)`$ at that density. If the original $`g_2(r)`$ is generated by a system with only pair interactions, then $`v(r;\rho )=w^{(2)}(r)`$ will be independent of density. If there are any higher order interactions in the original system which influence the structure, then this equivalence will only hold for $`lim_{\rho 0}v(r;\rho )=w^{(2)}(r)`$. At finite densities $`v(r;\rho )`$ must change since the structure is no longer equal to the one generated by $`w^{(2)}(r)`$ alone. Therefore $`v(r;\rho )`$ must be density dependent.

In fact, this is what we found in two previous papers, where we used the hypernetted-chain (HNC) approximation from liquid state theory to extract $`v(r;\rho )`$ from computer simulations of the $`g_2(r)`$’s between the CM’s of SAW polymers. For completeness, we show these effective pair interactions in Fig. (3). As expected, there is a clear density dependence. Without going into much detail about the inversion of $`v(r;\rho )`$ from $`g_2(r)`$, we do want to point out that the process can be very subtle. As illustrated in Fig. 4, the $`g_2(r)`$’s generated at $`\rho =\rho ^{}`$ by $`v(r;\rho =0)`$ and $`v(r;\rho =\rho ^{})`$ are very similar. Any technique to derive $`v(r;\rho )`$ from $`g_2(r)`$ must be significantly more accurate than the difference between the $`g_2(r)`$’s shown in the figure. The accuracy of the techniques we use has been discussed in Ref. , and will be analyzed in much more detail in another publication.

Any approximation which correctly reproduces the pair-correlations will also predict the correct thermodynamics through the compressibility equation. For our density dependent $`v(r;\rho )`$ this is indeed the case, since we found good agreement between the equation of state (EOS) $`\mathrm{\Pi }/\rho `$ generated by the effective potentials in Fig. 3 and the EOS of the underlying SAW polymer solution. In contrast, the $`v(r;\rho =0)`$ potential underestimates the EOS, and we find mean-field fluid behavior $`\mathrm{\Pi }/\rho \rho `$ at large $`\rho `$ instead of the correct $`\mathrm{\Pi }/\rho \rho ^{1.3}`$ scaling. So, even though the $`\rho =0`$ potential results in pair-correlations $`g_2(r)`$ that are similar to the true $`g_2(r)`$’s, the effective thermodynamics can differ significantly. The difference arises from the neglected many-body interactions, as discussed in the Appendix.

### B Understanding the density dependence of the effective pair-potential

Given the success of the density-dependent pair interaction in describing pair-correlations and thermodynamics, we next turn to the question of whether the density-dependence of $`v(r;\rho )`$ can be directly understood from the density independent many-body interactions.

Within the HNC approximation, the following expression due to Reatto and Tau and also Attard

$`v(r_{12};\rho )`$ $`=`$ $`w^{(2)}(r_{12})\rho {\displaystyle \left(e^{w^{(3)}(r_{12},r_{13},r_{23})}1\right)}`$ (14)
$`\times g_2(r_{13};\rho )g_2(r_{23};\rho )d𝐫_3,`$
describes the density dependence of the pair-potential that would reproduce the true pair-correlations induced by the 2-body and 3-body potentials. This is a generalization of earlier expressions and neglects terms of order $`𝒪(\rho ^2)`$ and higher. In the literature it has mainly been applied to the Axilrod-Teller interaction for rare-gas fluids, where it works remarkably well, see e.g. and references therein.

Figure 5 highlights the density dependence by plotting $`v_{ex}(r;\rho )=(v(r;\rho )w^{(2)}(r))/\rho `$. For clarity, we have replaced the rather noisy data by spline fits. For densities $`\rho /\rho ^{}<1`$ the curves are close to each other suggesting that the roughly linear density dependence in Eq. (14) holds true. For larger densities into the semi-dilute regime, $`v_{ex}(r;\rho )`$ becomes smaller in magnitude and the density dependence becomes non-linear. This non-linearity is not unexpected, since Eq. (14) neglects higher order terms in $`\rho `$, as well as the effects of 4-body and higher order interaction terms.

We can go even further and directly calculate the triplet induced density dependent term $`v_{ex}(r;\rho )`$ by substituting Eq (9) into Eq. (14) to obtain

$`v_{ex}(r_{12},\rho )`$ $`=`$ $`{\displaystyle \left(\underset{\rho 0}{lim}\frac{g_3(r_{12},r_{13},r_{23})}{g_2(r_{12})g_2(r_{13})g_2(r_{23})}1\right)}`$ (16)
$`\times g_2(r_{13};\rho )g_2(r_{23};\rho )d𝐫_3.`$
The $`g_2(r;\rho )`$, in contrast to the $`g_2(r)`$ in the first term in (16), are defined at the density of interest. Evaluating this integral is difficult, because the term between brackets can become very small. We use a direct MC procedure, where two polymer coils are held with their CM a distance $`r_{12}`$ apart while we integrate over the position of the third particle. In order to ensure that the integral converges it is crucial to use the $`g_2(r)`$ at $`\rho =0`$ (i.e. those between the brackets in Eq. 16), from the simulation itself, by calculating it on the fly. This is necessary to avoid small errors in the radial distribution function which build up during the integration over the volume. The $`g_2(r;\rho )`$ at finite density are known from previous calculations. The resulting $`v_{ex}(r;\rho )`$ for $`\rho =0`$ is plotted in Fig. 5. The results for finite density do not differ by very much.

In conclusion, the density dependence is mainly caused by 3-body interactions, at least in the dilute regime. At semi-dilute densities, higher order many-body interactions may come into play.

## V Many-body correlations

If $`g_2(r)`$ is generated by a system with only pair-potentials, then the exact inversion of $`g_2(r)`$ at any density will reproduce the exact pair potential. For such a system, the inverted pair potential can be used in principle to determine all higher order correlation functions of the original system. In other words, for systems with only pair-interactions, the pair-distribution function $`g_2(r)`$ contains enough information to generate all higher order correlation functions!

If a system has 3-body or higher order interactions, then our $`v(r;\rho )`$ still exactly reproduces the pair-correlations. But, as we shall demonstrate in this section, it can no longer exactly reproduce the higher order correlations. Nevertheless, we will show that the differences are not very large for the case of SAW polymers at the densities we study.

### A Bond angle distribution from 3-body correlations

Calculating and comparing the full 3-body correlation functions would be very cumbersome for many of the same reasons that it is difficult to map out the full 3-body interaction. Therefore, we resort to a reduced picture where a subset of the variables are integrated out . One popular measure of the 3-body interactions is the bond angle distribution function, defined as

$`b(\theta ,r_c)`$ $`=`$ $`8\pi ^2\rho ^2N{\displaystyle _0^{r_c}}{\displaystyle _0^{r_c}}g_3(r_{12},r_{13},(r_{12}^2+r_{13}^2`$ (18)
$`2r_{12}r_{13}cos\theta )^{1/2})r_{12}^2r_{13}^2sin\theta dr_{12}dr_{13},`$
where $`N`$ is a dimensionless normalization constant. This integral sums over all triplets within a cutoff radius $`r_c`$ from the central particle and determines the distribution of the angle $`\theta `$ in these triplets. We calculated the bond angle distribution for both the SAW simulations and the effective pair potentials for different cutoff radii $`r_c`$ as shown in Fig. 6. The effective potentials $`v(r;\rho )`$ reproduce this measure of the 3-body correlations remarkably well. Since for an ideal gas the bond angle distribution exactly follows a sine curve, dividing the bond angle distribution by $`\mathrm{sin}\theta `$ highlights the deviations from ideal behavior. In Fig. 7 we show the renormalized bond angle distribution. The differences between the curves are now clearer. The absolute deviations from the sine-like behavior are largest at small $`\theta `$ because the particles repel each other and triplets with small $`\theta `$ will be relatively rare. In the case of a hard sphere systems this correlation hole would be even more pronounced. It is remarkable how well the bond order distribution follows the ideal sine curve for the larger angles. At a cutoff radius of $`r_c=0.5`$ the distribution becomes flatter, reflecting the broad flat top of the repulsive Gaussian shaped pair potential.

Instead of determining the distributions from explicit simulations, we can also substitute the Kirkwood superposition approximation

$$g_3(r_{12},r_{13},r_{23})g_2(r_{12})g_2(r_{13})g_2(r_{23})$$
(19)
into Eq. 18 and calculate the integral directly by using the radial distribution functions from previous simulations. This approximation is also included in Fig. 6 and Fig. 7 and turns out to be very accurate, except for $`\theta 0`$ and $`\theta \pi `$ (see Fig. 7) where the simulations are prone to large statistical errors, due to the vanishing volume of the available phase-space.

The bond-angle distribution is not very sensitive to differences in the full 3-body correlations (see e.g. Ref. ), partially because it is an integrated quantity. An example of this is given in Fig. 8, where the bond angle distributions of the effective potentials $`v(r;\rho =0)`$ and $`v(r;\rho =\rho ^{})`$ are compared for the same density $`\rho =\rho ^{}`$. Clearly, there is hardly any difference between the distributions.

### B Deviations from Kirkwood superposition for 3-body correlations

A more sensitive measure of triplet correlations is the deviation from the Kirkwood superposition approximation, Eq. (19), which we define as:

$$G_3(r_{12},r_{13},r_{23})=\frac{g_3(r_{12},r_{13},r_{23})}{g_2(r_{12})g_2(r_{13})g_2(r_{23})}.$$
(20)
Since our $`v(r;\rho )`$ exactly reproduces the $`g_2(r)`$, this expression should highlight any differences between the true $`g_3`$ and the $`g_3`$ arising from our effective potential picture. To simplify, we limit ourselves for a given $`r_{12}`$ to triplet configurations for which $`r=r_{23}=r_{13}`$ (i.e. isoceles triangles).

First, we compare in Fig 9 the $`G_3(r)`$ at $`\rho =\rho ^{}`$ generated by $`v(r;\rho =0)`$ and by $`v(r;\rho =\rho ^{})`$. Just as we found for the bond angle distributions, the $`G_3(r)`$ are very similar even though the potentials are different. We already showed that the $`g_2(r)`$ are not very different either (see Fig. 4) , so that the same now holds for the full $`g_3(r_{12},r_{23},r_{13})`$.

Next, we turn to a comparison between the true $`G_3^{\mathrm{SAW}}(r)`$ derived from explicit simulations of our SAW polymer system and the $`G_3^{\mathrm{eff}}(r)`$ of the effective potentials at $`\rho =\rho ^{}`$ . As can be seen in Fig. 9, our effective pair-potential $`v(r;\rho )`$ does not exactly reproduce the SAW 3-body correlation function. The trends are similar, but the deviation from superposition of the SAW polymers, $`G_3^{\mathrm{SAW}}(r)`$, is consistently larger than the same quantity generated by the effective potentials $`G_3^{\mathrm{eff}}(r)`$, especially if $`r<R_g`$.

For systems with an explicit three-body interaction the Kirkwood superposition approximation is sometimes written as:

$`g_3(r_{12},r_{13},r_{23})`$ $``$ $`g_2(r_{12})g_2(r_{13})g_2(r_{23})`$ (22)
$`\times \mathrm{exp}\left[w^{(3)}(r_{12},r_{13},r_{23})\right].`$
This is exact in the $`\rho 0`$ limit, as can be seen from Eq. (9). Note that in this same limit the three-body correlations induced by the $`v(r;\rho )`$ reduce to the simpler Kirkwood superposition approximation $`g_3(r_{12},r_{13},r_{23})g_2(r_{12})g_2(r_{13})g_2(r_{23})`$, demonstrating explicitly that in contrast to the 2-body correlations, $`v(r;\rho )`$ cannot exactly reproduce the 3-body correlations if there is a 3-body interaction present! In fact, we have shown explicitly that two systems with identical pair correlations, namely our effective potential system and the original SAW system, can have differing triplet correlations.

These arguments also suggest that a simple approximation, namely $`G_3^{SAW}(r)G_3^{\mathrm{eff}}(r)\mathrm{exp}[w^{(3)}(r)]`$, can shed some light on the differences observed in Fig 9. Since $`w^{(3)}(r)`$ is negative for equilateral triangle configurations, as illustrated in Fig 2, it is perhaps not surprising that roughly speaking $`G_3^{SAW}(r)>G_3^{eff}(r)`$ for the isoceles triangle configurations plotted in Figs 9, at least in the region where $`w^{(3)}(r)`$ is non-zero. Unfortunately the statistical errors in this region are very large, making a quantitative comparison difficult, but the deviation is certainly of the same order as would be expected from an extra factor $`\mathrm{exp}[w^{(3)}(r)]`$ (compare with Fig.2).

Our very simple approximation illustrates how deviations from Kirkwood superposition originate from two effects: (1) Deviations induced by correlations generated by the pair-potentials alone. These have been studied in great detail for hard-sphere systems, see e.g.. (2) Deviations induced primarily by 3-body potentials. In practice, of course, these two effects are somewhat entwined, especially at higher densities. Nevertheless splitting the two effects can shed light on the origin of three-body correlations. In particular, it suggests that while an effective pair-potential $`v(r;\rho )`$ that exactly reproduces the $`g_2(r)`$ can partially reproduce deviations from superposition of type (1), it will fail for deviations of type (2).

Since the 3-body and higher order interactions between the CM of polymer solutions are not very strong, the total $`g_3(r)`$’s are still remarkably well reproduced by the $`v(r;\rho )`$, especially when integrated quantities such as the bond-angle correlations are considered. However, for systems where 3-body interactions are strong, such as liquid Si or liquid Ga, one cannot expect the same success from effective pair potentials. Very similar conclusions were stressed by Evans in the context of reverse Monte-Carlo simulations.

## VI conclusions

Integrating out the monomeric degrees of freedom to obtain a description based on effective potentials between polymer CM’s is a useful coarse-graining technique for polymer solutions. Because simulations can be performed to high accuracy, the lessons learned here should be applicable to a much broader range of coarse-graining schemes.

In particular, we showed that the free energy of the polymers can be expanded in a series of state independent many-body effective potentials. The terms in the series oscillate as $`(1)^n`$, and become smaller in absolute magnitude for increasing $`n`$. The scaling theory developed in the Appendix confirms these ideas, and can be used to extend them to arbitrary order $`n`$.

A parallel description of the coarse-grained polymer solution was developed in terms of an effective state (density) dependent pair potential $`v(r;\rho )`$, which exactly reproduces the pair-correlations and, in an average way, includes all the higher order terms in the many-body free energy expansion. The density-dependence of this effective pair-potential can be largely understood from the direct influence of the density-independent three-body interactions.

The three-body correlations are also well described by this effective pair-potential picture. If the bond angle distribution is used, the differences between the full SAW polymer system and our effective pair-potential picture are almost indistinguishable. However, this is not a good “order parameter” for measuring deviations in three-body distribution functions, since there is almost no difference between results from the full simulations and those produced with the Kirkwood superposition approximation. When we use a more direct measure of the deviations from Kirkwood superposition, small differences between the effective two-body and the full SAW triplet distributions can be measured. These arise primarily from the direct effect of the three-body potential $`w^{(3)}(r)`$, and illustrate a more general point, namely that an effective two-body interaction can never simultaneously reproduce the two and three-body correlations in a system with many-body interactions exactly. Since for polymers in a good solvent these many-body interactions are relatively weak, a coarse-grained description based on effective pair-interactions works remarkably well, at least for the dilute and the beginning of the semi-dilute regime. Whether this success can be extended deeper into the semi-dilute or into the melt-regime, or even for polymers in poor solvents, remains to be seen, and will be the subject of future investigations.

###### Acknowledgements.

AAL acknowledges support from the Isaac Newton Trust, Cambridge, PB acknowledges support from the EPSRC under grant number GR$`/`$M88839, We thank H. Löwen and C. von Ferber for useful discussions, and R. Finken and V. Krakoviak for a critical reading of the manuscript.

## A Relative strength of $`n`$ body terms from scaling theory

In a recent paper, von Ferber et al.  used scaling theory and simulations to calculate the triplet interaction for star-polymers and found an attractive interaction with a relative strength of about $`11\%`$, very similar to our 3-body results. The natural choice for the position coordinate of a star-polymer is not the CM, but its midpoint. However, we will see that for estimating relative contributions, the difference between our CM representation and the midpoint representation is not that important.

Here, we apply the star-polymer scaling theory to estimate the relative contributions of $`w^n(0)`$ to all orders in $`n`$. We specialize to linear polymers, which can be seen as star-polymers with only two arms. We first note that the partition function for $`n`$ polymers with their mid-points constrained to be at a distance $`rR_g`$ apart scales as:

$$Z_n(r)r^{\theta _n},$$
(A1)
in the limit $`r/R_g0`$. Here $`\theta _n`$ is the contact exponent which in turn can be written as

$$\theta _n=n\eta _2\eta _{2n}.$$
(A2)
where the $`\eta _f`$ are the scaling exponents for a star-polymer with $`f`$ arms. These are tabulated for two different renormalization group calculations in, and can also be approximated by a simpler expression

$$\eta _f0.335307f^{3/2}$$
(A3)
which is expected to become more accurate for larger $`f`$. Comparing the different approximations gives an indication of the accuracy of the scaling theory results.

The probability of finding $`n`$ polymers with their midpoints a distance $`rR_g`$ apart can be found from the partition functions since:

$$\underset{\rho 0}{lim}\mathrm{ln}g_n^{mid}(r)\mathrm{ln}\left[Z_n(r)\right]\theta _n\mathrm{ln}(r/R_g).$$
(A4)
Here $`g_n^{mid}(r)`$ is the $`n`$-body distribution function for the midpoints of $`n`$ polymers. These scaling theories can be developed for polymers in the midpoint representation in part because of the analogy of the overlap of $`n`$ polymers to a $`2n`$ arm star-polymer. Such an analogy is not as obvious for the CM representation, hampering the derivation of a CM based scaling theory. On the other hand, MC simulations for overlap in the CM representation are much easier than for the midpoint representation because the probability $`g_n(r)`$ is much larger than $`g_n^{mid}(r)`$ in the small $`r`$ limit, in fact $`g_n^{mid}(0)=0`$ while the CM $`g_n(0)`$ is finite. Gathering good statistics for the mid-point representation is thus much slower than for the CM representation. It nevertheless seems reasonable to expect that the relative strengths of different orders of the interactions in the CM representation are similar to those of the midpoint representation, suggesting that the relative probability of full overlap of the CM of $`n`$ polymers scales as:

$$\underset{\rho 0}{lim}\left(\frac{\mathrm{ln}g_n(0)}{\mathrm{ln}g_m(0)}\right)\underset{r0}{lim}\underset{\rho 0}{lim}\left(\frac{\mathrm{ln}g_n^{mid}(r)}{\mathrm{ln}g_m^{mid}(r)}\right)\frac{\theta _n}{\theta _m}.$$
(A5)
In Table I we confirm this ansatz by comparing $`\mathrm{ln}g_n(0)/\mathrm{ln}g_2(0)`$ calculated by direct simulations for $`n=3,4,5`$ with $`\theta _n/\theta _2`$ calculated by three different versions of the scaling theory. Clearly, all three scaling theories and the simulations agree reasonably well with each other, giving us confidence to proceed.

Armed with this approximate equivalence of the $`g_n`$’s in the two representations, we can now fruitfully compare the strength of the $`n`$ body interaction $`w^{(n)}(0)`$ for CM particles to an expression derived from the scaling theory, by using Eq. (A2) and Eq. (A4) to recursively rewrite Eq. (12) as:

$$w^{(n)}(0)\widehat{w}^{(n)}=(1)^n\underset{m=2}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{n}{m}\right)(1)^m(m\eta _2\eta _{2m}).$$
(A6)
Here $`\widehat{w}^{(n)}`$ is the coefficient of the full (divergent) midpoint-midpoint interaction $`w_{mid}^{(n)}(r)\widehat{w}^{(n)}\mathrm{ln}(r/Rg)`$ in the limit $`r0`$. The values of these terms relative to the $`n=2`$ term are compared in Table I with direct simulation results for $`n=3,4,5`$, and again good agreement is obtained. For the scaling theory each higher order term is opposite in sign and smaller than the previous one. Our simulations show the same trend for the sign, but the error bars on the $`n=5`$ simulations are still too large to ascertain that its magnitude is less than the $`n=4`$ term.

The strength of the $`n`$ body term contribution to the free energy relative to all lower order terms can also be found from the arguments above. In this case the pre-factors cancel and we find:

$$\mathrm{\Delta }w^{(n)}(0)\mathrm{\Delta }\widehat{w}^{(n)}=\frac{\widehat{w}^{(n)}}{\theta _n\widehat{w}^{(n)}}$$
(A7)
Once again the simulations and the scaling theory agree quite well both for the magnitude and the sign of the different terms.

## B Convergence of the full free energy series

The good agreement found between the scaling theory and the simulations for $`n=3,4,5`$ suggests that we can evaluate trends for arbitrary order $`n`$ from the scaling theory. We use the simple expression in Eq. (A3) which gives $`\eta _f`$ for arbitrary $`f`$, to simplify Eq. (A6) to

$$\widehat{w}^{(n)}=0.9484(1)^n\left(\underset{m=1}{\overset{n}{}}\left(\genfrac{}{}{0pt}{}{n}{m}\right)(1)^{m+1}m^{3/2}\right).$$
(B1)
The sum between brackets is positive for all $`n`$, and becomes smaller for each subsequent term. For example for $`n=3`$ it is $`0.29`$, while for $`n=10,000`$, it is has dropped to $`0.009`$. This implies that at all orders, the $`n`$-body interaction $`w^{(n)}(0)`$ oscillates in sign as $`(1)^n`$ and decreases in absolute magnitude with increasing $`n`$. The absolute magnitudes of the relative contributions $`\mathrm{\Delta }w^{(n)}(0)`$ decrease even faster with $`n`$.

We expect that the absolute magnitude of each $`n`$ body term is largest at full overlap. However, the fact that each higher order term decreases in absolute magnitude, or even that the terms oscillate in sign at $`r=0`$, does not necessarily imply that the full free energy expansion of Eq. (1) will converge quickly for all configurations $`\{𝐫_i\}`$ of the CM. For one thing, for an $`N`$ particle system, the number of $`n`$-tuple coordinates grows as the binomial $`N!/(Nn)!n!`$. Therefore, even though the magnitude of the $`n`$ body terms $`w^{(n)}(𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_n})`$ becomes smaller with increasing $`n`$, the number of such terms at each order $`n`$ in the free energy series Eq. (1) increases almost exponentially with increasing $`n`$. Of course for a given configuration $`\{𝐫_i\}`$, not only the magnitude, but also the range of each term $`w^{(n)}(𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_n})`$ is needed to decide how many, if not all $`N`$ terms, are needed to calculate the full free energy to a desired accuracy. As mentioned before, calculating the complete dependence of $`w^{(n)}(𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_n})`$ on all possible $`n`$-tuple CM coordinates $`\{𝐫_{i_1},𝐫_{i_2}\mathrm{}𝐫_{i_n}\}`$ is usually an impossible task for higher order $`n`$.

Nevertheless, some progress can be made by looking at the special case of a set $`\{𝐫_i\}`$ where all $`N`$ CM’s are at the same point. This is simply the situation studied in the previous section by scaling theory, for which the free energy is given by

$$(N,V,\{𝐫_i\})^{(0)}(N,V)=\mathrm{ln}g^N(0)\theta _N.$$
(B2)
Ignoring for the moment the volume term, who’s contribution is expected to be negligible we can ask the question: How well is the free energy described by taking only pair-interactions terms into account? Since there are $`\frac{1}{2}N(N1)`$ pair terms, the free energy with only pair-terms taken into account scales as $`N(N1)w^{(2)}(0)N^2`$ while the true free energy scales as $`\theta _NN^{3/2}`$. Including only the pair terms heavily overestimates the free energy, a result also found by von Ferber et al.. In other words, if one truncates the coarse-grained free energy series in Eq. (1) then certain configurations $`\{𝐫_i\}`$ will not be properly accounted for.

However, if one were to use this coarse-grained free energy series to study a given polymer solution, then the relative probability for finding a configuration at complete overlap of the CM would be very small. On the other hand, for any given order $`n`$, the probability of finding $`n`$-fold overlap increases with increasing density, so that one would expect to need more and more terms in the coarse-grained free energy series as the density increases.

Another way to measure the effects of truncating the coarse-grained free energy series in Eq. (1) would be to compare it to the excess free energy per particle for the original polymer solution. This can be found by averaging $`(N,V,\{𝐫_i\})`$ over all configurations $`\{𝐫_i\}`$ to obtain the total free energy $`F(N,V)`$, as done in Eq.(4). Including only pair interactions in our CM representation of polymer solutions casts the problem into that of finding the free energy of a mean-field fluid for which the excess free energy per polymer scales as:

$$\frac{F^{ex}}{N}\rho $$
(B3)
while the true excess free energy of a polymer solution in the semi-dilute regime follows from scaling theory

$$\frac{F^{ex}}{N}\rho ^{1/(3\nu 1)}\rho ^{1.3}.$$
(B4)
Whereas including only the pair interactions in an average over all configurations $`\{𝐫_i\}`$ in the semi-dilute regime underestimates the total free energy, taking into account only pair terms for single configuration $`\{𝐫_i\}`$ at complete overlap of the CM overestimates the coarse-grained free energy. This implies that the special configurations where the coarse-grained free energy series in Eq. (1) breaks down do not have a strong influence on the thermodynamics.


format: text

text:
# In search of an evolutionary coding style

## 1 Background

Let us as a background cite three different sources.

The first is J. Madox’s comment on page 376 in .

> “In genetics for example, the task of understanding the functions of all the 100,000 human genes will require a much greater effort than that involved in their identification, and by a factor 10 or more.”

Chris Adami from Caltech made, in his survey talk at Renaissance Technologies in Stony Brook 10/27/98 about artificial life, a brief remark about the quality of the evolved program codes in his avida set up. He said something like this:

> “The codes that are evolved will eventually be almost totally unreadable. Things are never used only once, but two or more times. It is a kind of a ‘madman’s’ code.”

In recent years, a lot of examples have been found where genes have been “reused” for different purposes during development. As an example, take the runt gene in Drosophila, which is used in sex determination, segmentation and central nervous system creation. In  the authors write:

> “As mature organisms we are composed of an astonishing array of diverse cell types—all derived from a single-celled zygote. When faced with the task of generating such cellular diversity in a reproducible fashion, how has the embryo chosen to respond? Recent work in a number of developmental systems has suggested that the embryo has employed two approaches. First, given finite resources, the embryo has efficiently chosen to reutilize a limited set of proteins in different temporal and spatial contexts to create cellular diversity. Second, the embryo has also chosen to install molecular redundancies to ensure the reproducibility of these patterns from individual to individual.”

How can we capture these comments about the style or quality of the computer code and the DNA, in a quantified manner? Can we do that in such a general manner that we will be able to use analogous quality measure both for carbon– and silicon based genetic codes?

### 1.1 Plan of the paper

We give a straightforward but general definition of a “code”, a general definition of the “style” of such codes using a given set of measures. We will also suggest such measures of characteristic features, such as some type of “madness”, the robustness, and the amount of reuse, of such codes. Finally, we will make some “baby–experiments” by running the avida program and analyze samples of the evolved code, generate two simulated programs in the same function class, and eventually stylistically compare the real evolved code with these simulated codes, all of this will be done in a C$`++`$ program. Some results will be graphically displayed at the end.

### 1.2 Future goals

One could then more systematically run the avida program (or something similar) and analyze the evolved codes stylistically with more realistic, and a higher number of comparasion codes. By changing parameters for the set up in avida, one might eventually capture some common features. By comparing the carbon based programming style and the silicon based style, it might be possible to find some common parts that would describe the natural programming style for evolution. That would in turn help us read carbon based code. And it would also give us some hints how to create more robust computer programs, by looking at how nature has solved such problems.

### 1.3 Questions

Is there an existing way to study the style or quality of the DNA? And is there any existing theory in the information sciences that deals with this on the silicon side? We are aiming at a complexity level higher than the usual information theory measures, such as the notion of entropy (see the comment in Section 6.2).

### 1.4 Acknowledgement

I would very much like to thank professor G. Thomsen for listen to my ideas and for his constant encouragement. I would also like to thank him for making me feel so much at home in his Xenopus lab, and teaching me the basics in frog handling and cultural behaviour in the life sciences.

I would also like to thank Dr D. Slice for his interest and all his help and suggestions.

B. Cohen has been a good source for me concerning the existing computer code complexity measures. He has also been a fruitful discussion partner.

I am grateful to S. Sutherland for reading an earlier version of this paper and giving me useful comments and pointing out errors and weaknesses.

Finally I would like to thank D. Brander for playing squash with me and for trying to teach me some grammar. (All the remaining language errors are naturally completely my own fault.)

## 2 A general code

We will give a general definition of a code as a string of generalized “letters” from a given alphabet, that when interpreted, will define a function. This interpretation is not unique, i.e. many different codes will produce the same function. That fact will give us a way to study classes of codes. That is, two codes are in the same class if their interpretation gives equivalent functions.

### 2.1 Codes

We will give a definition of a code by using an underlying alphabet, and a generalized interpretor.

###### Definition 2.1

Let us define a code to be a finite string of “letters” taken from an “alphabet”, $`𝒜`$,

$$\text{Code}=\{\alpha _i\}_{i=1}^k,\text{ where }\alpha _i𝒜,$$

such that when the code is interpreted, the code will represent a well defined function (or a process), $`\text{Code}_jf_j`$, with a domain $`D_{f_j}`$ such that for all inputs, $`xD_{f_j}`$, to the interpretation of the code, will give $`f_j(x)`$ as the output111Note that we here consider a function in a general sense, i.e. not necessarily numerical..

### 2.2 Classes of codes

From the above definition of codes, we see that different codes can have the same function representation, see Example 2.5 below.

Let us therefore introduce the following classification. Let the code $`\text{Code}_f`$ and the codes $`\text{Code}_{f_i}`$ have the functions $`f`$ and $`f_i`$ as their representations.

###### Definition 2.2

Let the code class with respect to the function $`f`$ be the following (infinite) set of codes.

$$𝒞_f=\{\text{Code}_{f_i}:f_i(x)=f(x),\text{ for all }xD_f\}.$$

###### Remark 2.3

If the interpretor is not able to produce a well defined function from the code Code, we say that Code is in the error class $`𝒞_\epsilon `$. That (huge) class can be viewed as the complement to all interpretable codes.

###### Question 2.4

With the suggested norm of $`𝒞_f`$ below, is $`𝒞_f`$ compact? Do there exist extremal codes?

###### Example 2.5

Let

$$f(x)=\frac{1}{x2},$$

$$g(x)=\{\begin{array}{cc}\frac{x+2}{x^24}\hfill & \text{ if }x2\hfill \\ \frac{1}{4}\hfill & \text{ if }x=2\hfill \end{array}$$

and

$$h(x)=\{\begin{array}{ccc}\hfill \frac{1}{4}& \text{ if }x=\hfill & \hfill 2\\ \hfill \frac{1}{2}& \text{ if }x=\hfill & \hfill 0\\ \hfill 1& \text{ if }x=\hfill & \hfill 1\\ \hfill 1& \text{ if }x=\hfill & \hfill 3\\ \hfill \frac{1}{2}& \text{ if }x=\hfill & \hfill 4\end{array}$$

Note that $`\text{Code}_g`$ and $`\text{Code}_f`$ are both in $`𝒞_h`$; and $`𝒞_f=𝒞_g𝒞_h`$.

###### Remark 2.6

One common way to view functions is as black boxes. Here we are interested of the internal structure of such black boxes performing equivalent tasks.

## 3 The style of a code

We will now give a description of a method to characterize different coding styles. Since different codes do not differ in function in the same class, we say that they differ in style. How can we characterize such style?

Given a set of measures on the codes, we propose a way to characterize style of a subset of codes in the code class as an extremal unit weight vector that will act as a stylistic “fingerprint”. We will also give an algorithm for “stylistic translations”, and an index which reveals how well a given subset represents a common style.

The methods we use apply very elementary mathematics, and maybe even more basic statistical methods. That will hopefully make it accessible to a wide scientific audience who are interested in “style”.

### 3.1 A measure on $`𝒞_f`$

Let us now study measures on $`𝒞_f`$. Let

$$\mu _i:𝒞_f[0,1].$$

Let us consider the following profile measure

$$𝝁=(\mu _1,\mu _2,\mathrm{},\mu _n)$$

of codes in $`𝒞_f`$. If a given measure has a range outside $`[0,1]`$, let us use the transformation

$$x\frac{x}{1+x}$$

to make it fit into $`[0,1]`$.

We define the following scalar measure.

$$\nu _𝐰(\text{Code}_g)=𝐰𝝁(\text{Code}_g),$$

where $`𝐰`$ is a normalized weight vector such that $`𝐰=1`$, for some norm $`||||`$. For example, let

$$𝐰=𝐰_p=\left(\underset{i=1}{\overset{n}{}}|w_i|^p\right)^{\frac{1}{p}},$$

where $`p1`$. As a default norm, let us use $`||||=||||_2`$.

### 3.2 Extremal weights

Let us use the above measure to try to capture a characterization of “style” of codes.

Suppose that we have two sample sets, $`A`$ and $`B`$, of codes in $`𝒞_f`$. We will try to find a style characterization of the codes in $`A`$ relative to $`B`$.

We can think of $`A`$ as the set of codes we are stylistically interested in and $`B`$ as a complementary environment.

Let us define a vector $`𝐮=(u_1,u_2,\mathrm{},u_n)`$ in the following way.

$$𝐮=\underset{a_iA}{}\underset{b_jB}{}(𝝁(a_i)𝝁(b_j)).$$
(1)
Let us now normalize $`u`$ to a unit vector.

$$𝐰^+=𝐰^+(A)=\frac{𝐮}{𝐮}.$$
(2)
Let us now study the random variable

$$X=\nu _𝐰(a_i)\nu _𝐰(b_j),$$
(3)
where $`a_i`$ is a randomly chosen code in $`A`$, with uniform probability222The probability to choose $`a_i`$ is $`1/\mathrm{\#}(A)`$., and $`b_j`$ is randomly chosen in $`B`$. Note that $`X`$ is dependent on the chosen $`𝐰`$.

###### Proposition 3.1

Let $`X`$ be the random variable defined in (3) and let $`𝐰^+`$ be the unit vector from (2). Then we have that picking $`𝐰=𝐰^+`$ will maximize the expected value of $`X`$, $`\text{E}(X)`$.

Proof: Let $`\mathrm{\#}()`$ denote the cardinality of a set and let $`M=\mathrm{\#}(A)\mathrm{\#}(B)`$. The expectation for a general weight vector $`𝐰`$ will be easily computed.

$$E(X)=\frac{1}{M}\underset{a_iA}{}\underset{b_jB}{}(\nu _𝐰(a_i)\nu _𝐰(b_j))=$$

$$=\frac{1}{M}\underset{a_iA}{}\underset{b_jB}{}(𝐰𝝁(a_i)𝐰𝝁(b_j))=\frac{𝐰𝐮}{M}.$$

Since $`𝐰=1`$ we have that

$$|E(X)|\frac{𝐮}{M}.$$

Let us now look at the special case $`𝐰=𝐰^+`$. We have then that

$$E(X)=\frac{𝐰^+𝐮}{M}=\frac{1}{M}\frac{𝐮}{𝐮}𝐮=$$

$$=\frac{𝐮^2}{M𝐮}=\frac{𝐮}{M}.$$

Thus we see that taking $`𝐰`$ to be $`𝐰^+`$ will maximize $`E(X)`$. $`\mathrm{}`$

### 3.3 A fingerprint on the hyper sphere

We can now view the unit vector $`𝐰^+(A)=𝐰^+`$ on the unit hyper333if $`n>3`$ sphere as a “fingerprint” of the style of codes in $`A`$ relative to $`B`$ with respect to the list of given measures in $`𝝁`$.

###### Remark 3.2

If there is a superset of codes $`B_1B`$ which is enough separated, then the fingerprint $`𝐰^+`$ of $`A`$ relative to $`B_1`$ can be expected to be a more refined characterization of the “style” in $`A`$ than the $`𝐰^+`$ of $`A`$ relative to $`B`$.

### 3.4 A universal character

One would hope that similar classes, $`𝒞_f`$, would give similar fingerprints $`𝐰^+`$ for different function codes created by the same code writing agent.

### 3.5 Stability with respect to the vector $`𝝁`$

There is a stability feature built into $`𝐰^+`$ in the sense that if you would like to find a stylistic quality in a group of codes, you try to find “relevant” measures in the vector $`𝝁`$. What happens if, in addition to your relevant measures, you also take a sequence of irrelevant measures (where, for example, the codes look more or less randomly distributed, or even just the same)? If the environment, i.e. $`B`$, is rich enough, then your profile will just be about zero at the tail, where the non-relevant measures are. That means you don’t have to be restrictive when you pick your measures — if some happen to be worthless, that will be taken care of by itself.

### 3.6 Is there a common style in $`A`$?

Given two subsets, $`A`$ and $`B`$ in $`𝒞_f`$, of codes, we have now a method for finding a common style in $`A`$, in relation to $`B`$, as a unit vector $`𝐰^+(A)`$. (As a special case, we can let $`B`$ be the complement of $`A`$ in $`𝒞_f`$. Then we can talk about the style of $`A`$.)

This process can be executed even if $`A`$ and $`B`$ are just randomly chosen subsets in $`C_f`$ where we can not expect to find any stylistic common features in $`A`$ in comparison to all the codes in $`B`$. How can we find out if $`A`$ really has a common style, in comparasion to $`B`$, that can be captured by the chosen measure profile $`𝝁`$?

Let us go back to Equation (3) and the random variable $`X`$. Let us also define a similar s.v. $`Y`$ such that

$$Y=\nu _𝐰(c_i)\nu _𝐰(c_j),$$

where $`c_i`$ and $`c_j`$ are randomly chosen, with equal probability, in $`AB`$. Note that the underlying $`𝐰`$ is $`𝐰^+(A)`$. From Proposition 3.1 above, the expected value of $`X`$ will then be maximal. Let us denote that value by $`m`$, i.e. let $`E(X)=m`$. From the proof of Proposition 3.1 we see that

$$m=\frac{𝐮}{\mathrm{\#}A\mathrm{\#}B}.$$

Hence, loosely speaking, $`m`$ is large if $`A`$ has a characteristic style that is captured by the measure profile $`𝝁`$. How large can $`m`$ get? Or in other words: how large can $`𝐮`$ get? Since all the measures $`\mu _k()`$ are bounded above by $`1`$ and below by $`0`$ we have that a component $`u_i`$ of the $`𝐮`$ vector is also bounded.

$$u_k=\underset{a_iA}{}\underset{b_jB}{}(\mu _k(a_i)\mu _k(b_j))$$

$$\underset{a_iA}{}\underset{b_jB}{}(10)=\mathrm{\#}A\mathrm{\#}B=M.$$

If we are using the standard norm $`||||=||||_2`$, we get that $`𝐮nM^2`$ and hence $`𝐮\sqrt{n}M`$.

Let us normalize $`m`$ to get the following index

$$\theta =\frac{m}{\sqrt{n}}.$$

We see that $`\theta [0,1]`$ and it will be closer to $`1`$ when the fingerprint is good.

On the other hand, due to the Remark 3.3 below, we might also need an invariant index.

Let us use the variances of the random variables $`X,Y`$ in the following way. Let

$$\sigma _A^2=E((Xm)^2),\text{ and let }\sigma _{AB}^2=E(Y^2).$$

Note that we immediately have that $`E(Y)=0`$. Let us suggest an index $`\eta `$ of how well the measure profile $`𝝁`$ can capture a common style, if it exists at all, in $`A`$ in comparison to the codes in $`B`$. Let

$$\eta =\eta (A,B,𝝁)=\frac{\sigma _{AB}^2}{\sigma _A^2}=\frac{E(Y^2)}{E(X^2)m^2}.$$
(4)
###### Remark 3.3

Note that if all the measures in $`𝛍`$ are multiplied by a factor $`k<1`$, i.e. $`\mu _ik\mu _i`$, then $`mkm`$, $`\sigma _A^2k^2\sigma _A^2`$, and $`\sigma _{AB}^2k^2\sigma _{AB}^2`$. Hence $`\eta `$ would not change, but $`\theta k\theta `$.

#### 3.6.1 Principal component analysis

A very illustrative, and very popular, technique when looking for connections in a multidimensional environment is the principal component analysis. It is a two dimensional diagram with coordinate axis the two eigenvectors with the largest eigenvalues of the covariance matrix $``$. Let $`C=AB`$ and let $`N=\mathrm{\#}C`$. We denote the codes in $`C`$ by $`c_j`$, then

$$_{i,j}=\text{Cov}((\mu _i(c_1),\mu _i(c_2),\mathrm{},\mu _i(c_N)),$$

$$(\mu _j(c_1),\mu _j(c_2),\mathrm{},\mu _j(c_N))).$$

For more details, see any book in multivariate analysis, for example .

If in such a diagram, the codes $`A`$ we are interested in are clearly separated from the $`B`$ codes, then we could say that there is a common style in $`A`$, and if that is the case and if the largest eigenvalue is considerably larger then the second one, then the fingerprint $`𝐰`$ and the first eigenvector should be close to each other. That is indeed the case in Figure 7.

#### 3.6.2 Cluster analysis

Another tool from the multivariate toolbox could be used to study the question about a common style in $`A`$. This method is based on an iteration fusion of close points until the desired number of subsets are obtained. To measure the closeness, we might pick our scalar measure $`\nu _{𝐰^+}`$.

In order to check if there is a common style in $`A`$, we can ask how well the points in $`A`$ are clustered.

### 3.7 Style translations by iterations

Using the above construction, let us here indicate an algorithm how to “stylistically translate” a code $`a`$ in $`A𝒞_f`$, into the style in $`B𝒞_f`$.

1. Compute

$$𝐯=\underset{b_iB}{}(𝝁(b_i)𝝁(a)).$$
2. Find the dominating component of the vector $`𝐯`$, i.e. let

$$v_m=\underset{1in}{\mathrm{max}}|v_i|.$$
3. Study the measure $`\mu _m`$ and stylistically rewrite $`a`$ such that $`\mu _m(a)`$ would increase approximately $`v_m`$ units (decrease if $`v_m<0`$).
4. Iterate the process until sufficient accuracy is attained. The accuracy is measured by $`𝐯`$.

Suppose the final accuracy is $`\delta `$ and denote the rewritten $`a`$ by $`a^{}`$. Let $`𝐰=𝐰^+(B)`$ and let $`Z`$ be the random variable

$$Z=\nu _𝐰(b_i)\nu _𝐰(a^{}),$$

for a randomly chosen $`b_i`$ in $`B`$. Now the expected value of $`Z`$ will be

$$E(Z)=\frac{_{iB}𝐰(𝝁(b_i)𝝁(a^{}))}{\mathrm{\#}(B)},$$

where $`\mathrm{\#}()`$ stands for the number of elements in the set. That is

$$E(Z)=\frac{𝐰𝐯}{\mathrm{\#}(B)}.$$

Now if we have the norm $`||||`$ as the default $`||||_2`$ we can use Cauchy–Schwarz’ inequality to get that

$$|E(Z)|\frac{\delta }{\mathrm{\#}(B)}.$$

In other words, the fingerprint $`𝐰^+`$ of $`B`$ would hardly “feel” the difference between $`a^{}`$ and the codes in $`B`$ if $`\delta `$ is small.

## 4 Different levels of Code

So far we have just been studying a code on a singular level. In this section we will describe a way to separate the codes into different levels. One could then ask questions about the styles on different levels. Are the styles similar even on different levels, etc?

We will view a code as a composition of lower level codes. Let us use the definition of emergence given in  on p. 518.

$$P\text{ is an emergent property of }S^2$$

$$$$

$$P\text{Obs}^2(S^2),\text{ but }P\text{Obs}^2(S_{i_1}^1)i_1.$$

Let us explain this more in detail, and similarly display a concrete example of a computer code, where the inputs are integers.

In this case, let $`S^0`$ be input values, e.g. $`S_i^0𝐍`$, and let $`\{\text{Int}^i\}_{i=0}^n`$ be a given sequence of sets of interactions, e.g. $`\text{Int}^0=\{+,\}`$, $`\text{Int}^1=\{,÷\}\text{Int}^0`$, $`\text{Int}^2=\{=,<,>\}\text{Int}^1`$, etc. Let $`\{\text{Obs}^i\}_{i=0}^n`$ be the related observational function, e.g. $`\text{Obs}^0(x)=`$ value of $`x`$ as an integer, $`\text{Obs}^2(x>y)=`$ true or false.

From the given sequences $`\{\text{Int}^i\}_{i=0}^n`$ and $`\{\text{Obs}^i\}_{i=0}^n`$ we get the higher order structures as a “reaction”, $`R`$, in the follwing manner.

$$S^1=R(S^0,\text{Obs}^0,\text{Int}^0),$$

$$S^2=R(S^1,\text{Obs}^1,\text{Int}^1),$$

see  for details.

We have then that in our example, $`1`$ is an emergent property of $`S^1=𝐙`$, $`1/2`$ is an emergent property of $`S^2=𝐐`$, and “true” is an emergent property of $`S^3`$, etc.

Given such sequence $`\{\text{Int}^i\}_{i=0}^n`$ we can define a code of degree $`k`$ as a consecutive string of the total code which has the property of $`\text{Obs}^k`$, i.e. if the code can produce an output that is an emergent property of $`S^k`$.

We can then view the final code of degree $`n`$, $`\text{Code}_n`$, as a composition of sub-codes of degree $`n1`$, $`\text{Code}_{n1}^i,`$ etc.

Note that Baas indicates this application in mentioning the word hyperalgorithms on p. 526 in .

How can one think of a good implementation of these functions? The choice of the interactions Int will give us a chance to find fine structures. What happens to our example if we start by $`\text{Int}^0=\{+\}`$ and then $`\text{Int}^1=\{\}\text{Int}^0`$ etc? Is there a “natural” choice of that interaction sequence for a given case?

## 5 Some measures

Let us list a couple of measures that would be useful to capture some of the features mentioned in the Section 1 above, and which also utilizes different levels of codes discussed above.

### 5.1 Spaghetti

Let us propose a kind of “code-madness-measure” using the above hierarchies of codes.

Let $`m_k`$ be the maximum numbers of $`\text{Code}_{k1}`$ codes in a $`\text{Code}_k`$. More precisely, Let

$$\eta _k^i=\mathrm{max}\{j:\text{Code}_{k1}^{i,j}\text{Code}_k^i\},$$

and let

$$m_k=\underset{i}{\mathrm{max}}\eta _k^i\text{ and }s_k=\underset{i}{}\eta _k^i.$$

Now, we define the “spaghetti length” at level $`k`$ to be

$$S_k(\text{Code}_n)=\frac{m_k}{s_k}\text{ and }$$

$$S(\text{Code}_n)=\underset{k}{\mathrm{max}}S_k(\text{Code}_n).$$

### 5.2 Reuse

Let $`0<kn`$

$$R_k^i(\text{Code}_n)=$$

$$\frac{\text{max }\mathrm{\#}\text{Code}_{k1}\text{Code}_k\text{ used }i\text{ times or more}}{s_k}.$$

Let us use the conventions $`R_k^2(\text{Code})=R_k(\text{Code})`$ and $`R_2^i(\text{Code})=R^i(\text{Code})`$.

### 5.3 Redundancy

A very important feature in evolutionary driven codes are their robustness. There are good measures of robustness given in the literature. One way to measure it is to to check the probability that the code “survives” a one point mutation.

The redundancy in a code of degree $`k`$, could be measured in the following way.

Let $`m`$ be the maximal number of subunit codes of type $`\text{Code}_{k1}`$ in the code of type $`\text{Code}_k`$, that can be taken away without affecting the output of the $`\text{Code}_k`$ code, and let $`n`$ be the total number of subunits. Let us then define the redundancy in code $`C`$ of type $`\text{Code}_k`$ by

$$Red(C)=\frac{m}{n}.$$

### 5.4 Brittleness

Another function that might be better to use in some ways would be a “brittleness” function that we define in the following way.

Let $`m`$ and $`n`$ be as above and let $`d(k)`$ be the number of $`\text{Code}_{k1}`$ codes in a $`\text{Code}_k`$ code, that when removed totally, or partially, destroys the code $`\text{Code}_k`$. The we define the brittleness of $`\text{Code}_k`$ as.

$$\text{Britt}(\text{Code}_k)=\frac{d(k)}{nm}.$$

That can be viewed in the following way. Let us think of $`\text{Code}_k`$ as a chain and the subunits as links in that chain. Then the length of the chain will be $`nm`$ and $`d(k)`$ will be the number of links that are can not be removed without breaking the chain.

Note also that in the special case when we either have one or two links for each step in the chain, that $`d(k)=n2m`$.

Note that the probability to survive a deletion of a randomly chosen subcode is $`1d(k)/n`$. That is, if there would just be one sublevel to the code, $`1d(k)/n`$ would be the usual robustness measure mentioned above.

## 6 Five applications

Since we are mainly interested in searching for a evolutionary universal coding style, we are of course interested in the two special cases when the code is a computer code, and when it is a sequence of the DNA. We have tried to make the above definitions of codes and style general enough to be able to deal with those cases.

As a by-product, we noted that this theoretical framework could also be applied to other areas where ‘style’ is essential. The applications that comes closes to our mind was art, music and literature. As a third application we will comment how the above stylistic fingerprint could be used for a common framework in “stylometry” in the study of authorship attribution. Our hope is that if a theory is not only applicable to the special cases it was aimed for under its construction, but also to a different case, it might be a sound approach in that theory.

A fourth possible application would be a stylistic investigation of the internal architecture of black boxes in the theory of neural nets, and its carbon based version—the brain.

The fifth application is about artificial life. We will make a small experiment in such an environment.

### 6.1 Computer code

By looking at the indicated toy-example above where we had natural numbers as inputs and $`+,`$ as the first interactions, it is not to hard to imagine that you would get the intuitive “usual meaning” increase in complexity in substructures in the code from lower order arithmetic operations, more complex functions, subroutines, program parts, the complete program.

Note that what usually is seen as a good programming style, that is separate codes into small, more or less, independent units, and not too many jumps back and forth, will give you a low $`S(\text{Code})`$; see Section 5.1.

Let us now mention three existing families of computer code measures. See also  and  for a detailed description of the two first examples and many others.

#### 6.1.1 Halstead’s Complexity Measures

In 1977 M. Halstead introduced a tool to measure the complexity of a computer program. It is perhaps the most well known measure of that kind.

The measure, or more precisely the family of five measures, is based directly on the code in the following way. Let $`n_1`$ be the number of distinct operators, $`n_2`$ the number of distinct operands, $`N_1`$ the total number of operators, and $`N_2`$ the total number of operands. From these numbers, the following measures are constructed.

| Program vocabulary | $`n=n_1+n_2`$. |
| --- | --- |
| Program length | $`N=N_1+N_2`$. |
| Difficulty | $`D=\frac{n_1N_2}{2n_2}`$. |
| Volume | $`V=N\mathrm{log}_2(n)`$. |
| Effort | $`E=DV`$. |

The Halstead’s measure seems to be a good candidate for measures in the stylistic search since it is only based on textual information and due to the fact that it has been used in many contexts and over such a long time and hence its properties are quite well known.

#### 6.1.2 McCabe’s Cyclomatic Complexity

In 1976, T. McCabe introduced a measure of the number of linearly independent paths through a computer program as a measure of the complexity of the code.

The cyclomatic complexity, $`CC`$, is computed in the following way. Let us study a schematic graph of the program and count the number of edges $`E`$, the number of nodes $`N`$, and the number of connected components $`c`$. Then

$$CC=EN+c.$$

If a program has a $`CC`$ higher than $`50`$ it is said to be unstable, since it is then “very likely” break down if it is altered.

The Cyclomatic Complexity is more of a measure of the inner logical complexity then the textual Halstead’s measure.

#### 6.1.3 GRASP

Let us describe an ongoing project at Auburn, Alabama, which address a numerical local measure of computer codes, and its display. See  for further information and down–loading of the program.

> “The overall goal of the GRASP project is to improve the comprehensibility of software. Thus, it is important to be able to identify complex areas of source code. The Complexity Profile Graph (CPG), a new graphical representation based on a composite of statement level complexity metrics provides the user with the capability to quickly recognize complex areas of source code. The CPG is significant in that it shows the complexity of a program unit as a profile of statement level complexity metrics rather than as a single, global metric.”

First the program code is parsed into non overlapping segments; then a series of measures, briefly described below, is applied to each segment, giving a local complexity measure of the program code.

The content complexity of a segment $`S`$ in the code is defined as

$$\eta (S)=log(\underset{TS}{}\text{Weight(T)}),$$

where $`T`$ are tokens in the segment $`S`$. For example, in , the weights for Ada 95 are given in the table below.

| Token Description | Symbol | Weight |
| --- | --- | --- |
| Logical operators | and, or, not, … | 1.5 |
| Comparison op. | $`<,>,=,<=`$, … | 1.5 |
| Left parenthesis | ( | 1.3 |
| Identifiers | var1, proc1, … | 1.0 |
| Others | $`+,,,/,),\mathrm{}`$ | 1.0 |

> “The context complexity provides a baseline level of complexity for segments of simple statements nested within a compound statement, which itself may be nested several levels deep. The complexity of a compound statement is based on three aspects: inherent complexity, reachability, and breadth.”

These three complexities are added, with weights, to obtain the context complexity.

Combining the content complexity and the context complexity, by a weighted sum, gives the profile metric $`𝝁(S)`$ for a segment $`S`$.

The local metric $`𝝁`$ can then be graphically presented as a histogram to give indications where the code is more complex. Thus that would help a programmer to point out code segments that would need some extra thought.

#### 6.1.4 Finding a coding style

Suppose that we would like to find a way to assign a specific style to a programmer. How could we do that?

Suppose he writes in C++. Then one might take a sample of his codes, look at his functions and subroutines. Go to other programmers and find as many functions as possible from that environment that one would like to be able to identify “our” programmer in the future.

Sort the codes into function classes $`𝒞_f`$ where our programmer has written $`a_f𝒞_f`$ and other programmers $`b_{f,j}B𝒞_f`$.

Let us now take a wide variety of measures of computer codes, such as those mentioned above, e.g. Halstead’s measures etc. Let us call the vector of those measures $`𝝁`$, where each $`\mu _i`$ is a function of a C++ code to the unit interval $`[0,1]`$.

Let use equation (1).

$$𝐮=\underset{f}{}\underset{b_{f,j}B}{}(𝝁(a_f)𝝁(b_{f,j})).$$

Now,

$$𝐰^+=\frac{𝐮}{𝐮}$$

will be the stylistic fingerprint of this programmer with respect to the vector of measures $`𝝁`$ and in the chosen environment. Note that $`𝐰^+`$ will be heavily dependent on the chosen reference codes, so it is important that that environment is chosen carefully.

### 6.2 DNA

It is hard to get a handle on the DNA code with usual computer code tools. As an example in  the authors study block entropies for DNA to find some indication of structure.

###### Definition 6.1

Let the length of the alphabet be $`\lambda `$. Then the (normalized) $`n`$–block entropy of a sequence is defined as

$$H_n=\underset{i=1}{\overset{\lambda ^n}{}}p_i^{(n)}\mathrm{log}_\lambda p_i^{(n)},$$

where $`p_i^{(n)}`$ is the probability of the $`i`$th combination of $`n`$ “letters”.

In the summary of  one can read the following.

> “Surprisingly, DNA sequences behave closer to completely random sequences than to written text. The very strict syntax of computer languages on the other hand is reflected by a very low average information content of its sub-strings.”

To us that means that the style of DNA sequences is not just close to chaos, but rather written in such a different style compared to written text, and even more different from computer codes.

The inputs, $`S^0`$ here are the amino acids, and one could take perhaps $`\text{Int}^0`$ to be “putting next to”, or sequencing.

The letters on the DNA level should be the four bases, the words are then three letters word coding for a amino acid. Add to the dictionary special start and stop words for the genes.

One could also study the situation on the amino acid level, i.e. one level up from the bare DNA code. That would give us 20 letters in an alphabet.444The alphabet is not unique, e.g. different spelling in example some bacterias and humans for some amino acids. The words would here be the genes.

#### 6.2.1 Hopes

One would hope that some hypothetical insight of a universal evolutionary driven style would give some hints how to better “read” the DNA code.

We would therefore try to find measures that are general enough in the below described computer experiment, so that some hypothesis about the DNA code could eventually be made.

There are indications that evolutionary driven code, may not be “optimal” in a basic sense, but might include peculiar turns and twists. See for example p. 180 in  where J. Madox discuss RNA Editing

> “The puzzle is to know why these changes, which are presumably advantageous to the organism, have not been incorporated in the gene themselves, thus avoiding the need for editing by way of afterthought—not to mention the need for a separate biochemical mechanism for carrying it out.”

He also addresses, on p. 203, the “junk code” in eukaryotic cells.

> “At the very least, this complication is an extra metabolic cost for eukaryotic cells. It is also potentially a source of error. What countervailing selective advantage can there possibly be in this arrangement?”

### 6.3 Literature

Let us now turn our attention to something different. Let us look at some examples in literature where “style” has been in the focus.

On p. 74 in  the author discuss computational differences between grammatical errors and stylistic weaknesses.

In , seven important problems with the existing authorship attribution studies are listed and discussed and some solutions proposed. The proposed solution to problem number three is to

> “study style in its totality. Approximately 1,000 style markers have already been isolated. We must strive to identify all of the markers that make up “style” — to map style the way biologists are mapping the genes.”

Furthermore, the suggested solution to problem five is to

> “Develop a complete and necessarily multi–faceted theoretical framework on which to hang all non-traditional authorship attribution studies.
>
> Publish the theories, discuss the theories, and put the theories to experimental tests.”

As an example of a suggested metric from the literature studies, one can take the Yule’s coefficient advocated among many others in . Let $`\{f_{ij}\}`$ represent the observed frequencies in a “two way contingency table” and let the Yule’s coefficient (see ) be defined as

$$Y=\frac{\sqrt{c1}}{\sqrt{c+1}},\text{ where }c=\frac{f_{11}f_{22}}{f_{12}f_{21}}.$$

#### 6.3.1 Stylometry: Finding the author

Let us now treat a hypothetical case, using the methods from Section 3, on the authorship attribution problem.

What do we mean by a code in this case? Let us look at Definition 2.1. We will interpret the “letters” as word taken from the “alphabet” which will in this case be a complete dictionary of the language in question. What about the functions?

Scarry gives in  a description how we can think about the reading process; see the following quote from the first chapter:

> “When we say ’Emily Brontë describes Catherine’s face,’ we might also say ’Brontë gives us a set of instructions for how to imagine or construct Catherine’s face.’ This reformulation is accurate if cumbersome, in that it shifts the site of mimesis from the object to the mental act.”

So, in this case we might think about the functions as descriptions, or simply constant functions, where the interpreter is the reader.555Here it is easy to see the crucial importance of the interpreter. I would for example be an extremely poor interpretor of a French text, even if I would get some picture of a story at the end. Examples of interpreted functions would be text where “boy meets girl”, or “prince meets ghost”. An even more refined version would be “boy meets girl described in an English sonnet”.

Suppose now that we want to test the hypothesis that author X has written a given sonnet $`s`$. Then we might gather all know sonnets written by X into the set $`A`$ as a subset of reference sonnets from that time period in $`S`$. Now $`a_iA`$ means sonnet number $`i`$ in $`A`$.

One wants to have as large environment $`S`$ as possible, but at the same time also as narrow as possible, i.e. from the same time period, etc. This has then to be chosen carefully and with great knowledge about the literature period and its authors. Let now $`B`$ be the set of comparison sonnets.

On the other hand when it comes to choose metrics, then we pick as many as it is numerically computable in reasonable time, which depends both on our time and our computer. For example, there has been some new interesting development using neural networks in authorship attributions, see for example . An even more exciting method was used in  where they used genetic algorithms in order to find the best features to measure. We can use such nets and genetically derived measures in our set of measures too!

Now, apply equation (1) to get the vector $`𝐮`$ which then is normalized to $`𝐰^+`$. That will be the stylistic finger print of author X, given the above constructed framework.

How good and reliable is this fingerprint? Suppose that $`A`$ is large and $`S`$ is rich, not only large but more or less complete with respect to the author representation. Then $`\eta `$ from equation (4) would be a good measure how reliable $`𝐰^+`$ is. We want $`\eta `$ to be large of course, but what is large enough?

To answer that we need to study the actual distributions of the sonnets evaluated by $`\nu _𝐰`$, make approximations and perform hypothesis tests.

The nice thing about this method is that one does not have to argue which measure should be applied. Just apply them all!

Talking about measures, there has been a very animated debate about the values of using computers in literary studies. R. Quiones said in  “Why don’t they simply read the plays?”. A literature expert reading a text is undoubtly very hard to beat when it comes to authorship attributions. But let us look at this situation having our suggested definition of measures of codes in mind. Isn’t the expert using a long array of measures, weighted together in an intricate and sometimes subconscious way? One measure could be: “X would never use that word in two consecutive sentences”. The measure would be the characteristic function of that event in the text, and the weight would be heavily negative. The more skill and familiarity the expert has about such tasks, the lager the array $`𝝁`$, and the more subtle the weighing process would be. It may seem like a trivialization to think in those terms, but since the task is very hard and complex, one would not expect that a computer would be programmed in the near future that would be generally better than a literature expert. Compare with the relatively straightforward problem to play chess.

### 6.4 Neural nets

P. Adams in  pointed out the possibility of strong evolutionary forces acting inside the brain in the process of learning. “Good” synapses will be rewarded by being strengthened, and “bad” will be punished by being weakened. What is good and bad are much more intricate and implicit qualities compared to the genetic evolution, but maybe evolutionary forces are in command in the learning process on a time scale of seconds instead of millenias. And if there is a universal coding style—that should then be found in the style of the programming of the neural nets too.

### 6.5 Artificial life

Avida, see  and , is a program in the Tierra class, see . Unlike Tierra it gives a natural 2 dimensional picture of an evolutionary process. Like Tierra, Avida is not a simulation of real carbon based DNA, but a real evolution in a silicon based (simple) world. We will now do an experiment in the Avida world.

## 7 A small experiment

Let us try some of the above concepts in a Avida experiment. We will deliberately try to make it as simple as possible in order to get some output in a straightforward way. But we will keep all the doors open to make variants and generalizations in the future.

We use Avida version 1.0.1 that is available on the web; see .

When running the program, you have the possibility to extract individuals and saving their data and their genetic code. To read more about this, see the documentation on  or even better—the book . The code and the data is saved to files such as 153-aagxs.

The data in 153-aagxs tells you, for example, that in addition to being able to replicate itself, the code also performs some other tasks. In this specific case, it takes input from a stack, performs a logical XOR twice, performs three NOTs, etc. For this, it is rewarded with a bigger time slice, and hence will reproduce (and survive) better.

### 7.1 The function class

Let us as the function $`f`$ defining the code universe $`𝒞_f`$, take a function that exactly performs the above described set of tasks, e.g. three NOTS etc.

As the subspace $`A`$ we simply take the generated avida code in 153-aagxs.

### 7.2 Comparison codes

To get a comparison environment, let us construct (or simulate) two, man–made codes that perform the operations in $`f`$.

We use a very simple approach and simulate two different codes in $`B`$. One with no loop, except for the self–copying loop, and one with as many loops as there are logical tasks to be done. Furthermore, we don’t actually write the codes but simulate the writing using a very rigid approach of making all the logical statements from combinations of NANDs666This is the base for the default reward list for completed tasks in time slicing; see the file task.set in the Avida package., and then in detail study the needed operations just in one NAND.

This approach will give us a list of operations, and operands needed to fulfill $`f`$. We can now try to look for style.

### 7.3 The measure

To start with something, we used the Halstead’s measures as our $`𝝁`$. That is, let

$$𝝁=(\text{vocabulary, length, difficulty, volume, effort}).$$

### 7.4 Simulations

The actual simulations is done by extracting creatures, i.e. individual programs, from a run of the Avida program, as described above. The extraction consists of some data, the table of performance, but also the program itself. We then use a $`C++`$ program to automatically analyze the extraction by first reading of the data and the performance table, and then find the Halstead’s measure vector $`𝝁`$ for the actual code. The $`C++`$ program then also uses the performance table to simulate the two different comparison codes described in Section 7.2 above and to calculate their Halstead’s measures respectively. The program then also computes $`𝐰`$ when $`A`$ has only the code of extracted individual, and $`B`$ consist of the two comparison codes. ( $`𝐰`$ for other combinations of $`A`$ and $`B`$ are also considered.) The results are then exported to a Maple file where one can more easily work with the output.

### 7.5 Preliminary results

Since $`\mathrm{\#}A`$ is one, and $`\mathrm{\#}B`$ is two, we can not talk about results in any statistical sense of course. Nevertheless, we can present some outcome that should be seen as a indicator of what eventually can be done. Even if we have a very small $`𝒞_f`$, since essentially every creature represents a unique list of performed tasks, which leads to essentially a unique function $`f`$; we have many classes $`𝒞_f`$, and we are free to change parameters, such as the initial random seed, the reward table for the tasks, etc. We can therefore compare profiles; see for example Figure 2.

Let us more in detail see what kind of output you can get by going back to our old friend, born after 19441 generations, 153-aagxs.

After feeding the file 153-aagxs into the C++ program, which disects the code, and simulates the two comparison codes, we get as outputs things like: The Halstead’s measures:

$$n_1=19,n_2=3,N_1=153,N_2=31$$

$$\text{difficulty}=98.1667$$

$$\text{volume}=820.535$$

$$\text{effort}=80549.2.$$

And we also get

$$\theta =0.0085549\text{ and }\eta =391865,$$

which does not really tell you much in this meager situation ($`\eta `$ is huge since the variance of $`X`$ is so extremely tiny for this special case). We can view $`𝐰^+`$ in a diagram; see Figure 3. (Remember that we normalized the measures into $`[0,1]`$ by the transformation $`\frac{x}{1+x}`$ in order to make the weights in $`𝐰`$ play in the same division.)