1 Introduction

In the strong-selection weak-mutation (SSWM) regime evolutionary dynamics reduces to an adaptive walk on what is known as a fitness landscape, the map from genotypes to fitness values (Fragata et al. 2019; de Visser and Krug 2014). For low mutation rates the nearly monomorphic population can be represented by a single majority genotype moving through the space of genotypes by individual mutations that fix with a probability depending on the fitness of the mutant relative to the parental genotype (Gillespie 1984; Orr 2002). Under strong selection, the movement of such a walker is additionally constrained towards increasing fitness values, making it an adaptive walk (Kauffman and Levin 1987). This limits the number of selectively accessible paths a population can take through the genotype space (Carneiro and Hartl 2010; Franke et al. 2011; Weinreich et al. 2005).

Here we investigate the impact that the mutational structure of the genotype space has on the number of evolutionary paths available to SSWM dynamics. We use a simple stochastic model for fitness landscapes known as the House-of-Cards (HoC) model, in which each genotype \({\varvec{g}}\) is assigned an i.i.d. continuous random fitness value \({F_{{\varvec{g}}}}\) (Kauffman and Levin 1987; Kingman 1978). Continuity of the distribution assures that ties in fitness values almost surely do not happen. Then a path is accessible if these fitness values are in increasing order. Accessibility therefore is a property purely of the ordering of the i.i.d. random variables. As a consequence accessibility is independent of the actual distribution chosen and we are free to choose any representative distribution. We will use the standard uniform distribution which has certain properties that make it easier to work with.

A genotype is made up of many individual sites or loci, which can be found in some given number of states called alleles and can be mutated individually. For simplicity we will assume that all loci have the same set of possible states. Therefore genotypes are sequences \({\varvec{g}}= {\left( { {g_{1}}, \ldots , {g_{L}} }\right) }\), with \(L\) determining the number of loci. Individual (point) mutations, which are the only ones to be considered here, mutate only one of the loci. The mutational structure of the system determines whether every state of one locus is able to mutate to any other or whether some restrictions apply. For example, whereas point mutations in the DNA sequence can mutate any nucleotide base into any other, the genetic code constrains the possible one-step transitions between amino acids. To accommodate general mutational structures we describe the loci by a simple directed graph.

Definition 1

An allele graph is a finite simple directed graph.

We will denote the allele graph under consideration \(\mathcal {A}\), its vertex set \(\mathbb {A}\), the size of its vertex set \(\mathfrak {A}\) and its adjacency matrix \(A\). The vertex set \(\mathbb {A}\) is the set of all alleles and the graph’s arrows indicate possible one-step mutations between alleles. The assumption of finiteness is not strictly necessary but allows for more focused proofs. Extension to infinite graphs is straight-forward with sufficient regularity properties (e.g. bounded degrees). Again for simplicity, we assume the allele graph to be the same on all loci with the vertices identified by natural numbers. We will give a justification for this restriction after introduction of the genotype space in the following paragraph.

Definition 2

A genotype space (over L loci with allele graph \(\mathcal {A}\)) is the Cartesian power graph \({\mathcal {A}^L}= \mathcal {A}^{\square L}\) (Martinsson 2018), where the Cartesian product \(\mathcal {U}\square \mathcal {V}\) of two directed graphs \(\mathcal {U}\) and \(\mathcal {V}\) is a graph over the Cartesian product of their vertex sets with an arrow from (uv) to \((u',v')\) iff either there is an arrow from u to \(u'\) in \(\mathcal {U}\) or there is an arrow from v to \(v'\) in \(\mathcal {V}\), but not both.

With this definition the genotype space’s vertex set forms the genotypes, with arrows between genotypes that can be reached via one-step mutations (Fig. 1). The fitness landscape constrains which of these arrows may be taken by an adaptive walker and we call the directed sub-graph of the genotype graph obtained by removing arrows which do not point towards increasing fitness the fitness graph (Crona et al. 2013) (Fig. 2). For conciseness we usually do not specify the number of loci and the allele graph for the genotype space under consideration, assuming these parameters to be \(L\) and \(\mathcal {A}\) instead.

Fig. 1
figure 1

Example of a genotype space as the Cartesian graph product of two allele graphs. The two allele graphs are shown on top with the genotype space below. While the second factor graph represents a locus with three possible alleles, all of which may mutate freely from one to another, the first factor graph represents a locus on which not all mutations between the alleles are considered possible. Specifically mutations between 0 and 2 must take an intermediate mutation through 1 and additionally the mutation from 0 to 1 is considered irreversible and does not allow backstepping to 0. Although in this work we define the genotype graph as the direct Cartesian power of a single allele graph, different allele graphs as shown here can still be modeled without loss of generality by assuming that \(\mathcal {A}\) is the disjoint sum of the individual graphs. Since the individual constituent graphs are not connected in this sum, this does not increase the accessibility (see Remark 1)

Fig. 2
figure 2

Example of a fitness graph generated from the genotype space in Fig. 1 according to the HoC model. The opacity of nodes indicates the randomly chosen fitness value. Only arrows representing mutations that were originally allowed in the genotype graph and also point towards increasing fitness remain, resulting in an acyclic directed graph on genotypes. The global minimum and maximum in this realization are (10) and (12) respectively and the latter is accessible from the former by multiple accessible paths, for example the direct one \((10) \rightarrow (12)\), but also \((10) \rightarrow (20) \rightarrow (22) \rightarrow (12)\). As a counter-example to accessibility consider (01) and (12). Although (12) has higher fitness than (01), it is not accessible from (01)

Remark 1

Although we assume the allele graph to be the same on all loci, this still implicitly covers situations in which different loci use different numbers of alleles, as long as the number of alleles is bounded by a constant. This can be seen by considering an allele graph constructed as the disjoint union of all allele graphs with size at most the constant bound. The resulting genotype space then also separates into many disjoint components, one of which is the original genotype space with the varying sequence of allele graphs. Because mutations between the components are impossible, we can consider, without loss of generality, this larger genotype space without affecting the accessibility property.

Definition 3

A fitness landscape on a genotype space \({\mathcal {A}^L}\) is an assignment of real numbers \({F_{{\varvec{g}}}}\), called fitness values, to each genotype \({\varvec{g}}\). The fitness graph of the fitness landscape is the directed graph over all genotypes with an arrow from \({\varvec{g}}\) to \({\varvec{h}}\) iff there is an arrow from \({\varvec{g}}\) to \({\varvec{h}}\) in \({\mathcal {A}^L}\) and \({F_{{\varvec{h}}}} > {F_{{\varvec{g}}}}\).

As seen in this definition we will use bold face for genotypes, i.e. vertices of the genotype space, while using normal face for alleles, i.e. vertices on the allele graph. Throughout, quantities defined over the allele graph will be written in normal face and analogous quantities defined over the genotype space will be written in bold face.

Definition 4

A House-of-Cards (HoC) model over a genotype space \({\mathcal {A}^L}\) is a random distribution over the set of all fitness landscapes on \({\mathcal {A}^L}\), such that each fitness value is chosen i.i.d. from a standard uniform distribution.

Definition 5

Given genotypes \({\varvec{g}}\) and \({\varvec{h}}\) on a fitness landscape over a genotype space \({\mathcal {A}^L}\),

  1. 1.

    a walk on \({\mathcal {A}^L}\) from \({\varvec{g}}\) to \({\varvec{h}}\) is called accessible if it is also a walk on the fitness graph of the fitness landscape and

  2. 2.

    \({\varvec{h}}\) is said to be accessible from \({\varvec{g}}\) if at least one such walk exists.

We write \({\varvec{Z}}_{{{\varvec{g}}}{{\varvec{h}}}}\) for the number of accessible walks from \({\varvec{g}}\) to \({\varvec{h}}\), considering it as a random variable over the HoC distribution.

Definition 6

The (HoC) accessibility of a genotype \({\varvec{h}}\) from a genotype \({\varvec{g}}\) on a genotype space is the probability that \({\varvec{h}}\) is accessible from \({\varvec{g}}\) in a HoC model over the genotype space.

In other words accessibility of \({\varvec{h}}\) from \({\varvec{g}}\) is the quantity \({\mathbb {P}\left[ {\varvec{Z}}_{{{\varvec{g}}}{{\varvec{h}}}} \ge 1\right] }\).

Our goal is to determine the accessibility between pairs of genotypes \({{\varvec{a}}_L}\) and \({{\varvec{b}}_L}\) defined on genotype spaces with \(L\) loci as \(L\) becomes large (with fixed allele graph). This question is in particular non-trivial if the directed distance \({{\varvec{d}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}}\) from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\) on the genotype space is of linear order in \(L\). Here we understand directed distance \({{\varvec{d}}_{{{\varvec{g}}}{{\varvec{h}}}}}\) between two genotypes \({\varvec{g}}\) and \({\varvec{h}}\) as the length of the shortest (not necessarily accessible) path from \({\varvec{g}}\) to \({\varvec{h}}\) on the genotype space. Of special importance for this question is the value \(\beta \), which we define as the fitness difference \({F_{{{\varvec{b}}_L}}}-{F_{{{\varvec{a}}_L}}}\).

Definition 7

For values \(0 \le \beta \le 1\), the \(\beta \)-(HoC-)accessibility of a genotype \({\varvec{h}}\) from a genotype \({\varvec{g}}\) on a genotype space is the probability that \({\varvec{h}}\) is accessible from \({\varvec{g}}\) in the HoC model, given that the HoC distribution is conditioned on \({F_{{\varvec{h}}}}-{F_{{\varvec{g}}}} = \beta \).

We typically still write \({\mathbb {P}\left[ {\varvec{Z}}_{{{\varvec{h}}}{{\varvec{g}}}} \ge 1\right] }\) to refer to \(\beta \)-accessibility. If there is ambiguity between the two in a given context, we add notation to indicate the conditioning.

It is known from previous work on the case of two alleles, \(\mathbb {A}= \{0,1\}\), and linear distance \({{\varvec{d}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}} \sim \delta L\), that there is a critical value \(\beta _\text {c}\), depending on \(\delta \), such that for constant choices of \(\beta \) above or below \(\beta _\text {c}\), asymptotically the probability of \({{\varvec{b}}_L}\) being \(\beta \)-accessible from \({{\varvec{a}}_L}\) converges to 1 or 0, respectively, as \(L\rightarrow \infty \) (Berestycki et al. 2016, 2017; Hegarty and Martinsson 2014; Li 2018; Martinsson 2015). The transition occurring at \(\beta = \beta _\text {c}\) has been referred to as accessibility percolation (Krug 2021; Nowak and Krug 2013). Apart from a computational study (Zagorski et al. 2016), so far accessibility percolation has been studied only for the biallelic case for which the genotype space \({\mathcal {A}^L}\) is the \(L\)-dimensional (binary) hypercube.

Remark 2

Results related to those presented here have been obtained in the context of first-passage percolation (Kistler and Schertzer 2020; Martinsson 2016, 2018). The accessibility percolation problem (with any continuous distribution) can be mapped to an equivalent first-passage percolation problem with uniformly distributed weights as described in Martinsson (2015). This mapping does however map accessibility percolation with fitness values on vertices to first-passage percolation with weights on vertices as well, while traditionally weights are put on edges in that context. Martinsson (2018) considers a first-passage percolation model that would map to the HoC accessibility percolation problem if fitness values were assigned to edges rather than vertices. We will adapt and extend his methods to directly resolve the specific accessibility problem introduced here without requiring the mapping to first-passage percolation.

Remark 3

Another way of looking at the HoC \(\beta \)-accessibility problem is to consider it as a Bernoulli percolation on a certain ensemble of orientations of \({\mathcal {A}^L}\). Alternatively to the definition of \(\beta \)-accessibility above, the \(\beta \)-conditioning in the HoC model can also be applied by conditioning on \({F_{{{\varvec{a}}_L}}} = 0\) and \({F_{{{\varvec{b}}_L}}} = \beta \). This follows immediately from the fitness values being i.i.d. uniform random variables. If all fitness values are increased by a constant value, taking their remainder modulo 1, the resulting distribution is unchanged, so that the \(\beta \)-accessibility must be unaffected by further conditioning on the initial fitness value. With this then only genotypes with fitness values below \(\beta \) are relevant in determining whether an accessible path exists. Furthermore, after removal of the ineligible genotypes, accessibility will depend only on the order of the fitness values on the remaining vertices. Because the fitness values are chosen i.i.d. this implies that \(\beta \)-accessibility in the original problem is equivalent to 1-accessibility after additional removal of each vertex (that is not \({{\varvec{a}}_L}\) or \({{\varvec{b}}_L}\)) with probability \(1-\beta \), i.e. a Bernoulli site percolation with rate \(\beta \). This perspective may be more suitable if one is to consider the effect of decreasing \(\beta \).

2 Results

In this section we present our general results for (mostly) arbitrary allele graphs \(\mathcal {A}\) of which specific applications will be demonstrated in Sect. 3. First we define some attributes of the problem more carefully and state the limitations imposed by our subsequent proofs.

2.1 Prerequisites

Our intent is to describe the limiting behavior of accessibility as \(L\rightarrow \infty \) for fixed allele graphs between genotypes in distances which are large, i.e. linear in \(L\). To make this setup precise we need to introduce some more notation. We also introduce some restrictions on the choices of the allele graph \(\mathcal {A}\) and the choices for the sequence of endpoints \({{\varvec{a}}_L}\) and \({{\varvec{b}}_L}\) in order to avoid pathological and non-converging behavior. These restrictions and assumptions are summarized in Definition 10.

First we note that since allele graphs were defined to be finite, there is a maximal degree \(\varDelta < \infty \) among all vertices. When extending the results to infinite allele graphs it should be required that such an upper bound on the degree still exists in order for most results to remain valid. If the degrees are not sufficiently bounded in an infinite allele graph, then the number of walks from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\) may become so large that the problem results in trivial accessibility.

As \(L\), and with it the graph \({\mathcal {A}^L}\), changes we need to define a sequence of endpoint pairs \(({{\varvec{a}}_L},{{\varvec{b}}_L})\) in \(L\). Correspondingly we also intend \(\beta \) to vary with \(L\) as we consider \(\beta \)-accessibility. For convenience the dependence of \(\beta \) on \(L\) is taken to be implicit and usually not reflected in notation. For general sequences of endpoints calculations may become tediously complex and so we impose a few restrictions described in the following on the sequences for which our results will apply. Because the position of a locus in the sequence of the genotypes does not matter, all relevant properties of the pair \(({{\varvec{a}}_L},{{\varvec{b}}_L})\) for a given \(L\) can be expressed as an integer-valued matrix:

Definition 8

The allele counting matrix of a pair of genotypes \(({{\varvec{a}}_L},{{\varvec{b}}_L})\) over a genotype space is the \((\mathfrak {A},\mathfrak {A})\)-matrix with entries

$$\begin{aligned} M_{{v_{}}{w_{}}} = \left| \{{a_{l}}= {v_{}}\wedge {b_{l}}= {w_{}}\;|\; l=1\cdots L\}\right| \end{aligned}$$
(1)

where \({v_{}},{w_{}}\in \mathbb {A}\) are alleles.

This matrix counts for each pair of alleles the number of loci on which the path is required to move from \({v_{}}\) to \({w_{}}\), thereby dividing out the permutation-symmetry of loci. A sequence of such matrices M in \(L\) is equivalent to a sequence of pairs \(({{\varvec{a}}_L},{{\varvec{b}}_L})\) up to the irrelevant symmetry of \({\mathcal {A}^L}\).

For each value of \(L\) we have then a matrix M corresponding to the pair \(({{\varvec{a}}_L},{{\varvec{b}}_L})\) of endpoints. Because we are interested in the behavior as \(L\) diverges, we want to focus on cases where this sequence of matrices is sufficiently well-behaved. Borrowing the terminology from phylogeny studies (Lockhart et al. 1994), we introduce the following object:

Definition 9

The divergence matrix for a sequence of pairs of genotypes \((({{\varvec{a}}_L},{{\varvec{b}}_L}))_{L}\) with allele counting matrices \(M^{(L)}\) on a corresponding sequence of genotype spaces with \(L\) alleles and fixed allele graph is the \((\mathfrak {A},\mathfrak {A})\)-matrix with elements

$$\begin{aligned} p_{{{v_{}}}{{w_{}}}} = \lim _{L\rightarrow \infty }\frac{M^{(L)}_{{v_{}}{w_{}}}}{L} \end{aligned}$$
(2)

for all \({v_{}},{w_{}}\in \mathbb {A}\), assuming all limits exist. We also define a sequence of divergence remainder matrices with elements

$$\begin{aligned} R^{(L)}_{{v_{}}{w_{}}} = M^{(L)}_{{v_{}}{w_{}}} - Lp_{{{v_{}}}{{w_{}}}}. \end{aligned}$$
(3)

We want this divergence matrix to exist and the remainder terms to be sufficiently well-behaved that we can make convergence statements about the accessibility question. Specifically we require the following:

Definition 10

A (well-behaved) accessibility setup is a sequence \(({\mathcal {A}^L})_L\) of genotype spaces with \(L\) loci on the same allele graph \(\mathcal {A}\) together with a sequence of pairs of genotypes \((({{\varvec{a}}_L},{{\varvec{b}}_L}))_L\) on these genotype spaces such that

  1. 1.

    for all \(L\): \(M^{(L)}_{vw} = 0\) whenever there is no non-zero length walk from v to w in \(\mathcal {A}\),

  2. 2.

    for all \(L\): \({{\varvec{a}}_L}\ne {{\varvec{b}}_L}\),

  3. 3.

    the divergence matrix p of \((({{\varvec{a}}_L},{{\varvec{b}}_L}))_L\) exists,

  4. 4.

    at least one off-diagonal element of p is non-zero.

The first requirement assures that there is always at least one candidate path from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\), so that the trivial case does not need to be considered and that it is actually possible to move on all loci. If there are only zero-length walks on a locus (i.e. start and end point are the same without any possibility of mutation), then effectively the locus cannot affect accessibility in any way, resulting in a pathological effective reduction of \(L\).

The second requirement assures that we do not need to consider the pathological case in which the \(\beta \)-conditioning cannot be satisfied.

The third requirement assures that the “direction” between the pair of endpoints is asymptotically well-behaved without oscillations which would need to be reflected in the critical \(\beta \).

The last requirement assures that the distance between the endpoints grows linearly in \(L\), which is the only limit we consider here.

In the following we are always working with one such implied accessibility setup. In order to succinctly state our results, we introduce the following quantity for pairs of alleles \({v_{}},{w_{}}\in \mathbb {A}\) and \(t> 0\):

$$\begin{aligned} \varGamma _{{{v_{}}}{{w_{}}}}{\left( {{t}}\right) }&= \ln {\left( {{\left( \textrm{e}^{{t}A}\right) _{{{v_{}}}{{w_{}}}}}}\right) }, \end{aligned}$$
(4)

with \(A\) the adjacency matrix of the allele graph. The exponential is a matrix exponential from which the element representing alleles \(({v_{}},{w_{}})\) is extracted, rather than the exponential of the element of the matrix. In addition to the quantity \(\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{t}}\right) }\), also its first two derivatives \(\varGamma '_{{{v_{}}}{{w_{}}}}{\left( {{t}}\right) }\) and \(\varGamma ''_{{{v_{}}}{{w_{}}}}{\left( {{t}}\right) }\) with respect to \(t\) will be important. We indicate the derivative with respect to the \(t\) argument by backticks as shown.

Proposition 1

\(\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{t}}\right) }\) on \(t\ge 0\) has the following properties:

  1. 1.

    If \({v_{}}\ne {w_{}}\) and there exists no walk from \({v_{}}\) to \({w_{}}\) on the allele graph, it is nowhere defined (or \(-\infty \) everywhere),

  2. 2.

    if \({v_{}}= {w_{}}\) and there exists no walk of non-zero length from \({v_{}}\) to \({w_{}}\) on the allele graph, it is 0 everywhere,

  3. 3.

    otherwise it is strictly monotonic increasing, converging to positive infinity as \(t\rightarrow \infty \) and to negative infinity as \(t\rightarrow 0\).

In case 1. we formally interpret \(\textrm{e}^{\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{t}}\right) }}\) as \({\left( \textrm{e}^{{t}A}\right) _{{{v_{}}}{{w_{}}}}}\), which would be zero everywhere.

Proof

These properties are direct consequences of the behavior of the matrix exponential. \(\square \)

We then define the same quantity for pairs of genotypes \({\varvec{v}},{\varvec{w}}\in {\mathbb {A}^L}\) as averages over the per-locus quantity:

$$\begin{aligned} {\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }&= {\left\langle {\varGamma _{{{v_{l}}}{{w_{l}}}}{\left( {{t}}\right) }}\right\rangle }_{l} \end{aligned}$$
(5)

where we use the general notation

$$\begin{aligned} {\left\langle {X_l}\right\rangle }_{l} = \frac{1}{L}\sum _{l=1}^LX_l \end{aligned}$$
(6)

to mean an average over l of the term \(X_l\) containing l as a variable. As we already did for genotypes, quantities acting on \({\mathcal {A}^L}\) will be written in boldface, while equivalent quantities acting on a single \(\mathcal {A}\) copy will be denoted in normal font-face. The canonical connection between the former and latter is averaging over loci. In particular, \({\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }\) can be interpreted as the exponential rate with which the expected number of accessible paths grows in \(L\), or more precisely the number of quasi-accessible walks, a concept which will be introduced in Sect. 4.

Proposition 2

In a (well-behaved) accessibility setup \({\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }\) satisfies the following properties on \(t> 0\):

  1. 1.

    \({\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }\) is strictly monotonic increasing in \(t\), diverging to \(\infty \) as \(t\rightarrow \infty \) and to \(-\infty \) as \(t\rightarrow 0\),

  2. 2.

    \({\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }\) is bounded from above by a continuous strictly monotonic increasing function of \(t\),

  3. 3.

    \({\varvec{\varGamma }}'_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }\) is bounded from below and above by continuous strictly positive functions of \(t\).

  4. 4.

    All higher derivatives of \({\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }\) are bounded from below and above by continuous functions.

Proof

The conditions on the allele counting matrix in an accessibility setup guarantee that all contributions to the mean in Eq. (5) fall under the last point of Proposition 1. Furthermore because the allele graph is finite and fixed, there is only a finite number of pairs \(({v_{}}, {w_{}})\) of loci. The properties in this proposition then follow immediately from \({\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }\) being a point-wise average of a finite number of continuous strictly monotic increasing funtions with the same properties and only the relative weighting dependent on \(L\). \(\square \)

As a consequence of the first property there exists exactly one value for each \(L\) at which \({\varvec{\varGamma }}_{{{v_{}}}{{w_{}}}}{\left( {{t}}\right) }\) becomes 0, which we denote \({\hat{\beta }}\). This is going to be our candidate for the threshold of \(\beta \)-accessibility.

We also require the following function with domain \(0 \le r,s\le 1\), which is a slight generalization of a function introduced by Martinsson (2018):

Definition 11

Martinsson’s function of an accessibility setup is the function on \([0,1]^3\) defined by

$$\begin{aligned} {\mathfrak {M}({s},{r},{\beta })} = {\left\langle {{\left\langle { \varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) } }\right\rangle }^{s,r}_{{x_{l}},{y_{l}}}}\right\rangle }_{l}, \end{aligned}$$
(7)

where \({{\bar{r}}}= 1-r\) and \({{\bar{s}}}= 1-s\) and \({\left\langle {\cdot }\right\rangle }^{s,r}_{{x_{l}},{y_{l}}}\) is the mean over \({x_{l}},{y_{l}}\in \mathbb {A}\) weighted by

$$\begin{aligned} \textrm{e}^{\varGamma _{{{a_{l}}}{{x_{l}}}}{\left( {{\beta {{\bar{s}}}r}}\right) } + \varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) } + \varGamma _{{{y_{l}}}{{b_{l}}}}{\left( {{\beta {{\bar{s}}}{{\bar{r}}}}}\right) }}. \end{aligned}$$
(8)

Explicitly written out it reads

$$\begin{aligned} {\mathfrak {M}({s},{r},{\beta })}&={\left\langle { \frac{ \sum _{{x_{l}},{y_{l}}\in \mathbb {A}}\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) }\textrm{e}^{\varGamma _{{{a_{l}}}{{x_{l}}}}{\left( {{\beta {{\bar{s}}}r}}\right) } + \varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) } + \varGamma _{{{y_{l}}}{{b_{l}}}}{\left( {{\beta {{\bar{s}}}{{\bar{r}}}}}\right) }} }{ \sum _{{x_{l}},{y_{l}}\in \mathbb {A}}\textrm{e}^{\varGamma _{{{a_{l}}}{{x_{l}}}}{\left( {{\beta {{\bar{s}}}r}}\right) } + \varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) } + \varGamma _{{{y_{l}}}{{b_{l}}}}{\left( {{\beta {{\bar{s}}}{{\bar{r}}}}}\right) }} }}\right\rangle }_{l}. \end{aligned}$$
(9)

For the case that \({y_{l}}\) is not reachable from \({x_{l}}\) or that \(s= 0\), the formula yields negative infinities for \(\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) }\). We assume that in this case formally the natural choice

$$\begin{aligned} \varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) }\textrm{e}^{\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) }} = 0 \end{aligned}$$
(10)

holds. Martinsson’s function can be interpreted as a refined version of \({\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{\beta }}\right) }\) in which walks are segmented into three parts, each of which traverses fitness spans of \(\beta {{\bar{s}}}r\), \(\beta s\) and \(\beta {{\bar{s}}}{{\bar{r}}}\) respectively, weighting the expected number of accessible walks on the middle segment by the number of accessible walks by which it can be reached from \({{\varvec{a}}_L}\) and \({{\varvec{b}}_L}\). In this way a positive value of Martinsson’s function relative to \({\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\) shows that the number of accessible walks is clustered around few initial and end segments with many alternative accessible middle segments, which in turn implies that the overall expected value is not a good indicator for the existence of at least one accessible walk. How this function arises from consideration of a variation of the second moment method is described in Sect. 6.

Proposition 3

\({\mathfrak {M}({s},{r},{\beta })}\) satisfies the following properties:

  1. 1.

    \({\mathfrak {M}({0},{r},{\beta })} = 0\),

  2. 2.

    \({\mathfrak {M}({1},{r},{\beta })} = {\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\),

  3. 3.

    \({\mathfrak {M}({1},{r},{{\hat{\beta }}})} = 0\).

Proof

The first two statements can be derived immediately by application of matrix multiplication. The third statement is an immediate consequence of the second given the definition of \({\hat{\beta }}\). \(\square \)

The objects introduced so far are dependent on \(L\) implicitly through the averaging process over loci. In order to be able to make statements about the limiting behavior, it is useful to consider the limits of these quantities as \(L\rightarrow \infty \). We use the non-\(L\)-dependent mean

$$\begin{aligned} {\left\langle {X_{{v_{}}{w_{}}}}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}} = \sum _{{v_{}},{w_{}}\in \mathbb {A}} p_{{{v_{}}}{{w_{}}}} X_{{v_{}}{w_{}}} \end{aligned}$$
(11)

as a replacement for \({\left\langle {X_{{v_{}}{w_{}}}}\right\rangle }_{l}\). As the number of alleles is finite, \({\left\langle {X_{{v_{}}{w_{}}}}\right\rangle }_{l}\) will converge to \({\left\langle {X_{{v_{}}{w_{}}}}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}}\) as \(L\rightarrow \infty \).

In particular we write \(\beta ^*\) for the limit of \({\hat{\beta }}\) as \(L\rightarrow \infty \). Similarly we write for the limit of Martinsson’s function

$$\begin{aligned} {\mathfrak {M}^*({s},{r},{\beta })}&= {\left\langle {{\left\langle { \varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) } }\right\rangle }^{s,r}_{{x_{l}},{y_{l}}}}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}}. \end{aligned}$$
(12)

Proposition 4

In a (well-behaved) accessibility setup \(\beta ^*= \lim _{L\rightarrow \infty } {\hat{\beta }}\) exists, is positive and satisfies \({\left\langle {\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{\beta ^*}}\right) }}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}} = 0\).

Proof

Per the requirements on the divergence matrix for an accessibility setup the weights in the averaging over pairs of loci in \({\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\) converge. Because the average is only over a finite number of terms, it consequently also converges. Continuity of the matrix exponential and the properties from Proposition 2 then immediately prove this proposition. \(\square \)

With the necessary quantities defined we can classify accessibility setups as follows:

Definition 12

If not only at \(s= 0\) and \(s= 1\), but everywhere in its domain \({\mathfrak {M}^*({s},{r},{\beta ^*})} \le 0\), then we say that the accessibility setup is of semi-regular type. Otherwise we say that it is of irregular type. If \({\mathfrak {M}^*({s},{r},{\beta ^*})} < 0\) holds strictly everywhere except at \(s= 0\) and \(s= 1\) and if additionally the derivative \(\partial _s{\mathfrak {M}^*({s},{r},{\beta ^*})}\) is not zero at \(s= 1\), then we say that the accessibility setup is of regular type.

For our statements and in the following proofs we make use of Landau notation with the usual meanings of \(\mathcal {O}{\left( {{\cdot }}\right) }\), \(o{\left( {\cdot }\right) }\), \(\omega {\left( {\cdot }\right) }\) and \(\varTheta {\left( {\cdot }\right) }\). In our notation of arithmetic terms and equations these symbols are stand-ins for some function in the respective class. The limit variable to which these symbols apply should be evident from context, but is usually \(L\rightarrow \infty \). Functions in these classes are not required to be non-negative. In particular e.g. \(|\varTheta {\left( {1}\right) }|\) is used to enforce positiveness of a term that is of constant (non-zero) asymptotic order in \(L\). If not stated otherwise, these symbols are assumed to be uniform in the sense that the functions represented depend only on the limit variable and model parameters, but not on other local variables.

2.2 Statements

We state our results in terms of (weak) threshold functions defined as follows.

Definition 13

Given a function \(f(L)\) we say that a sequence \(c_L\) in \(L\) is a \(f(L)\)-threshold function for \(\beta \)-accessibility if for all \(g(L) = \omega (f(L))\)

  1. 1.

    \((c_L+|g(L)|)\)-accessibility has a non-zero limit inferior as \(L\rightarrow \infty \) and

  2. 2.

    \((c_L-|g(L)|)\)-accessibility has a zero limit superior as \(L\rightarrow \infty \).

We refer to the first condition as the upper side and the second condition as the lower side of the threshold. In other words, \(c_L\) determines the asymptotic transition from zero accessibility to non-zero accessibility if we condition the fitness difference between initial and final genotype, with a window of uncertainty of the same order as \(f(L)\). In particular if \(c_L\) is a \(f(L)\)-threshold for some \(f(L) = o(1)\), then the limit of \(c_L\) is the critical value \(\beta _\text {c}\).

Remark 4

The notion of threshold chosen here is weak in the sense that it doesn’t imply a transition from zero to one, but only from zero to some non-zero probability. We do not think that our results are actually restricted to this weak bound and we expect that arguments analogous to those made in Martinsson (2018) may be used to extend our weak threshold result to a strong threshold with \(\liminf _{L\rightarrow \infty }{\mathbb {P}\left[ {\varvec{Z}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} > 0\right] } = 1\), but we did not pursue this improvement here.

We can now state our main theorems:

Theorem 1

In an accessibility setup with the common notation used in previous definitions and the sequence

$$\begin{aligned} c_L= {\hat{\beta }}- {{\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) }}^{-1}\frac{\ln L}{L}, \end{aligned}$$
(13)
  1. 1.

    \(c_L\) satisfies the lower side condition for a \(\frac{1}{L}\)-threshold function and

  2. 2.

    \(c_L\) also satisfies the upper side condition and therefore is a \(\frac{1}{L}\)-threshold function if the accessibility setup is of regular type.

The second part of this theorem in particular implies that \(\beta _\text {c}= \beta ^*\) for the regular type, meaning that \(\beta ^*\) is indeed the critical value.

Theorem 2

In an accessibility setup of irregular type \(\beta _\text {c}> \beta ^*\) if \(\beta _\text {c}\) exists.

In other words, the irregular type does not have \(c_L\) as given in Theorem 1 as a threshold function.

Lastly we consider the length of accessible paths at the critical point:

Theorem 3

In an accessibility setup at the candidate threshold function from Theorem 1, i.e. with \(\beta = c_L + \frac{\eta }{L}\) for some constant \(\eta \), the probability that all \(\beta \)-accessible paths have length in the interval \({\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) } {\hat{\beta }}L\pm g(L)\sqrt{L}\) converges to 1 for every function \(g(L) = \omega (1)\).

In other words, if there are accessible paths at the candidate threshold function, then they have up to fluctuations of order \(\sqrt{L}\) length \({\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) }{\hat{\beta }}L\) and since Theorem 1 guarantees that the candidate is actually the threshold function for regular setups, this then implies that the critical paths in the regular setup are of that length as well.

The lower side of the threshold functions in Theorem 1 can be derived directly from a consideration of the expected number of (quasi-)accessible walks and an application of Markov’s inequality. This approach will be explained in Sect. 4, where we also introduce the notion of quasi-accessibility as a tool to simplify the counting of accessible paths. In addition, the first-moment approach allows us to prove Theorem 3 by consideration of the expected values separated by walk length (Sect. 4.3).

To prove the upper side of the threshold function, it is necessary to bound a higher moment of the expected number of (quasi-)accessible walks in relation to the mean. In particular, using a generalized version of the second moment method, it is sufficient to bound moments of the form

$$\begin{aligned} {\mathbb {E}\left[ \frac{{\varvec{Z}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}}{{\mathbb {E}\left[ {\varvec{Z}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }}\ln \frac{{\varvec{Z}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}}{{\mathbb {E}\left[ {\varvec{Z}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }}\right] } \end{aligned}$$

to show asymptotic boundedness of the accessibility away from zero. The evaluation of this expected value will follow the general ideas used by Martinsson (2018) to bound for every given (quasi-)accessible focal walk the number of other (quasi-)accessible walks, through the deviating arcs on the focal walk that generate all such other walks. In the mean taken over \({x_{l}}\) and \({y_{l}}\) in Martinsson’s function (9), the focal walk is represented by the walk sequence

$$\begin{aligned} {a_{l}}\rightarrow {x_{l}}\rightarrow {y_{l}}\rightarrow {b_{l}}\end{aligned}$$
(14)

and the corresponding three \(\varGamma \)-terms in the weights, while the deviating arcs are represented by the additional term corresponding to \({x_{l}}\rightarrow {y_{l}}\) over which the average is performed. Because Martinsson considers a model that corresponds to putting weights on edges rather than nodes, our calculations need to be adjusted accordingly (see Sect. 6).

The lower bound on \(\beta _\text {c}\) for the irregular case in Theorem 2 is again obtained following an approach used by Martinsson, by considering walks through pairs of edges \(({\varvec{x}},{\varvec{x}}')\) and \(({\varvec{y}},{\varvec{y}}')\), applying Markov’s inequality separately, and union bounding the resulting probability to improve on Markov’s inequality from the total expected number of (quasi-)accessible walks (see Sect. 5).

2.3 Asymptotic form

The theorems as stated in the previous section are dependent on \({\hat{\beta }}\) and \({\left\langle {\cdot }\right\rangle }_{l}\) averages, which are \(L\)-dependent quantities. From the definition of an accessibility setup we do however know that \({\hat{\beta }}\) converges to \(\beta ^*\) and that averages of the form \({\left\langle {\cdot }\right\rangle }_{l}\) are asymptotically of the form \({\left\langle {\cdot }\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}}\), both of which are \(L\)-independent quantities. Depending on the specific choice of sequence of pairs \((({{\varvec{a}}_L},{{\varvec{b}}_L}))_L\), the rate of convergence for these quantities may however differ and add additional significant terms in the threshold function, which we detail in this section.

Proposition 5

In Theorem 1 the sequence \(c_L\) can be replaced by some sequence

$$\begin{aligned} c_L= \beta ^*- {{{\varvec{\varGamma }}^{\prime *}}}^{-1}\frac{\ln L+ (1+o(1))\sum _{{v_{}},{w_{}}\in \mathbb {A}}R_{{v_{}}{w_{}}}\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) }}{L} \end{aligned}$$
(15)

for a suitable choice of o(1), and under this change in Theorem 3 there is a suitable choice of o(1) so that the interval becomes

$$\begin{aligned} {{\varvec{\varGamma }}^{\prime *}}\beta ^*L+ (1+o(1))\sum _{{v_{}}{w_{}}} R_{{v_{}}{w_{}}} {\left( {\varGamma '_{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) } - \frac{{{\varvec{\varGamma }}^{\prime \prime *}}}{{{\varvec{\varGamma }}^{\prime *}}}\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) }}\right) } \pm g(L)\sqrt{L} \end{aligned}$$
(16)

with

$$\begin{aligned} {{\varvec{\varGamma }}^{\prime *}}= {\left\langle {\varGamma '_{{{a_{l}}}{{b_{l}}}}{\left( {{\beta ^*}}\right) }}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}}, {{\varvec{\varGamma }}^{\prime \prime *}}= {\left\langle {\varGamma ''_{{{a_{l}}}{{b_{l}}}}{\left( {{\beta ^*}}\right) }}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}}. \end{aligned}$$
(17)

Proof

The leading order of the distance between the endpoints as \(L\rightarrow \infty \) is given by the sum of off-diagonal terms of \(p_{{{v_{}}}{{w_{}}}}\):

$$\begin{aligned} {{\varvec{d}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}} = \delta L+ o{\left( {L}\right) } \end{aligned}$$
(18)

where \(\delta = \sum _{{v_{}}\ne {w_{}}}p_{{{v_{}}}{{w_{}}}}\). The definition of an accessibility setup enforces that \(\delta >0\), so the value of \(\beta ^*\) will be positive, i.e. not zero, and then we can expand \({\hat{\beta }}\) around \(\beta ^*\) in \(L\), using that all \(\varGamma \)-terms are bounded by constants from below and above per Propositions 1 and 2:

$$\begin{aligned} 0 = {\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) }&=\mathrel {}{\left\langle {\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{\beta ^*}}\right) }}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}} + \frac{1}{L}\sum _{{v_{}},{w_{}}\in \mathbb {A}}R_{{v_{}}{w_{}}}\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) } \end{aligned}$$
(19)
$$\begin{aligned}&\quad +\,{\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta ^*}}\right) }{\left( {{\hat{\beta }}-\beta ^*}\right) } + \mathcal {O}{\left( {{{\left( {{\hat{\beta }}-\beta ^*}\right) }^2}}\right) }. \end{aligned}$$
(20)

The term \({\left\langle {\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{\beta ^*}}\right) }}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}}\) is zero by definition of \(\beta ^*\) and \({\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta ^*}}\right) }\) is bounded from below by a positive constant, so solving for \({\hat{\beta }}\) results in

$$\begin{aligned} {\hat{\beta }}&= \beta ^*- (1+o(1))\frac{1}{{\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta ^*}}\right) }L}\sum _{{v_{}},{w_{}}\in \mathbb {A}}R_{{v_{}}{w_{}}}\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) } \end{aligned}$$
(21)
$$\begin{aligned}&= \beta ^*- (1+o(1))\frac{{{\varvec{\varGamma }}^{\prime *}}^{-1}}{L}\sum _{{v_{}},{w_{}}\in \mathbb {A}}R_{{v_{}}{w_{}}}\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) }. \end{aligned}$$
(22)

Inserting this into the candidate threshold function (13) we obtain the alternative threshold function in the proposition.

Similarly, expanding \({\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{c_L + \frac{\eta }{L}}}\right) }\) in Theorem 3 gives, up to terms of order \(\frac{\ln L}{L}\) or smaller

$$\begin{aligned}&{{\varvec{\varGamma }}^{\prime *}}+ {{\varvec{\varGamma }}^{\prime \prime *}}{\left( {c_L - \beta ^*}\right) } + \frac{1}{L}\sum _{{v_{}},{w_{}}\in \mathbb {A}}R_{{v_{}}{w_{}}}\varGamma '_{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) } + \cdots \end{aligned}$$
(23)
$$\begin{aligned}&\quad ={{\varvec{\varGamma }}^{\prime *}}+ {\left( {1+o(1)}\right) }\frac{1}{L}\sum _{{v_{}},{w_{}}\in \mathbb {A}}R_{{v_{}}{w_{}}}{\left( {\varGamma '_{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) } + \frac{{{\varvec{\varGamma }}^{\prime \prime *}}}{{{\varvec{\varGamma }}^{\prime *}}}\varGamma _{{{v_{}}}{{w_{}}}}{\left( {{\beta ^*}}\right) }}\right) } + \cdots \end{aligned}$$
(24)

The remaining terms can contribute at most \(\ln L\) to the walk length, which would be subsumed by the interval length of order \(\sqrt{L}\). \(\square \)

In general, the candidate critical value does not depend on the non-linear corrections in the behavior of \(({{\varvec{a}}_L},{{\varvec{b}}_L})\), but the leading correction to the critical value changes if the non-linear corrections are of order \(\ln L\) or higher.

Remark 5

If \(\sum \limits _{{v_{}},{w_{}}\in \mathbb {A}}|R_{{v_{}}{w_{}}}| = \mathcal {O}{\left( {{1}}\right) }\), then the sequence \(c_L\) in Proposition 5 reduces to

$$\begin{aligned} c_L= \beta ^*- {{{\varvec{\varGamma }}^{\prime *}}}^{-1}\frac{\ln L+ \mathcal {O}{\left( {{1}}\right) }}{L}. \end{aligned}$$
(25)

This condition describes the situation in which \({{\varvec{a}}_L}\) and \({{\varvec{b}}_L}\) are, up to discretization error, separated in a well-defined linear “direction” of a fixed allele counting matrix as \(L\) increases. In this case the \(\mathcal {O}{\left( {{1}}\right) }\) contribution is irrelevant since \(c_L\) represents a \(\frac{1}{L}\)-threshold function, so that additional contributions of order \(\frac{1}{L}\) lie within the threshold window. In this linear separation case the threshold function is described fully by the two quantities \(\beta ^*\) and \({{\varvec{\varGamma }}^{\prime *}}\), both of which were derived from the averages over the matrix exponential of the allele graph, weighted by the divergence matrix.

3 Applications

3.1 Complete graph

The simplest application is to the complete graph on \(\mathfrak {A}= |\mathcal {A}|\) alleles, as seen in the top-left of Fig. 3, which leads to genotype spaces known as Hamming graphs. By symmetry, in this case there are only two choices for the initial and final allele on a locus, either \({a_{l}}= {b_{l}}\) or \({a_{l}}\ne {b_{l}}\). Therefore the accessibility setup can be fully determined by just the relative distance \(\delta \), which is then also the relative Hamming distance. As shown in Martinsson (2018), this accessibility setup is (for converging \(\delta > 0\)) always of regular type for the complete graph. One obtains

$$\begin{aligned} {\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }&= -\ln \mathfrak {A}- \beta + \delta \ln {\left( {-1 + \textrm{e}^{\mathfrak {A}\beta }}\right) } + {\bar{\delta }}\ln {\left( {\mathfrak {A}-1 + \textrm{e}^{\mathfrak {A}\beta }}\right) }, \end{aligned}$$
(26)

where \({\bar{\delta }} = 1-\delta \). In the biallelic case \(\mathfrak {A}= 2\) the condition \({\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) } = 0\) reduces to the relation \(\sinh ({\hat{\beta }})^\delta \cosh ({\hat{\beta }})^{\bar{\delta }} = 1\) which was first conjectured in Berestycki et al. (2016) and proved in Li (2018) and Martinsson (2015). At full distance \(\delta = 1\) without any variation of \(\delta \) with \(L\), \({\hat{\beta }}=\beta ^*\) and

$$\begin{aligned} {\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }&= -\ln \mathfrak {A}- \beta + \ln {\left( {-1 + \textrm{e}^{\mathfrak {A}\beta }}\right) }. \end{aligned}$$
(27)

The values of \(\beta ^*\) and \({{\varvec{\varGamma }}^{\prime *}}\) for small \(\mathfrak {A}\) are shown in Table 1 and numerically obtained values for \(\beta ^*\) as a function of distance are presented in Fig. 4.

Fig. 3
figure 3

Allele graph structures described in Sect. 3. Top left: Complete graph on four alleles with backmutations to the wild-type allele 0. Top right: Complete graph on four alleles with backmutations to the wild-type allele 0 removed. Bottom: Path graph on four alleles. In each case the intended initial (wild-type) and final alleles as used in this section are indicated by the labels \({a_{l}}\) and \({b_{l}}\)

Table 1 Results for the complete allele graph with 2–4 loci, as well as 21 loci, at full distance \(\delta = 1\)
Fig. 4
figure 4

\(\beta ^*\) as a function of the relative distance \(\delta = \lim _{L\rightarrow \infty }\frac{{{\varvec{d}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}}}{L}\) for the complete allele graph with \(2-6\) alleles. On the complete graph the distance specifies the choice of \({{\varvec{a}}_L}\) and \({{\varvec{b}}_L}\) up to irrelevant symmetries, so that this figure captures the full behavior

In general \(\textrm{e}^{\beta ^*}\) is the unique positive solution of the polynomial equation

$$\begin{aligned} {\left( {\textrm{e}^{\beta ^*}}\right) }^{\mathfrak {A}} - \mathfrak {A}\textrm{e}^{\beta ^*} - 1 = 0. \end{aligned}$$
(28)

For \(\mathfrak {A}\ge 5\) the solution of this equation cannot be expressed in closed form, however it can be expanded around \(\mathfrak {A}\rightarrow \infty \) as

$$\begin{aligned} \beta ^*&= \frac{\ln \mathfrak {A}}{\mathfrak {A}} + \frac{1+\ln \mathfrak {A}}{\mathfrak {A}^2} + \mathcal {O}{\left( {{\frac{\ln \mathfrak {A}}{\mathfrak {A}^3}}}\right) }, \end{aligned}$$
(29)
$$\begin{aligned} {{\varvec{\varGamma }}^{\prime *}}&= \mathfrak {A}+ \mathcal {O}{\left( {{\frac{1}{\mathfrak {A}}}}\right) }, \end{aligned}$$
(30)
$$\begin{aligned} \beta ^*{{\varvec{\varGamma }}^{\prime *}}&= \ln \mathfrak {A}+ \frac{1+\ln \mathfrak {A}}{\mathfrak {A}} + \mathcal {O}{\left( {{\frac{\ln \mathfrak {A}}{\mathfrak {A}^2}}}\right) }. \end{aligned}$$
(31)

As the number of alleles increases, accessibility increases and the required fitness difference between the start and end point decreases. In fact this quantity vanishes to zero for \(\mathfrak {A}\rightarrow \infty \). At the same time the length of accessible walks close to the critical fitness difference increases, but slowly. The minimal length of a path covering the full distance \({{\varvec{d}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}}\) is \(L\), and hence \(\beta ^*{{\varvec{\varGamma }}^{\prime *}}-1\) is the fraction of mutational reversions (where a mutated locus reverts to the allele it carried in the initial genotype \({a_{l}}\)) and sideways steps (where a mutation occurs to an allele that is part of neither the initial nor the target genotype) (Wu et al. 2016). The fraction of all alleles on a given locus that appear along an accessible path close to the critical point is given by \(\mathfrak {A}^{-1} \beta ^*{{\varvec{\varGamma }}^{\prime *}}\) which decreases with increasing \(\mathfrak {A}\). Zargorski et al. carried out simulations of this model, giving \(\beta ^*\) with two digit precision for different values of \(\mathfrak {A}\) (Zagorski et al. 2016). Their results match the values derived here up to \(\pm 0.01\).

3.2 Complete graph without return to the wild type allele

We can modify the complete graph slightly to disallow mutations back to the allele that was present in the initial genotype (the wild type allele), while still allowing mutations between all other allele, see top-right of Fig. 3. In this case the expressions simplify significantly to

$$\begin{aligned} {\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }&= \ln \frac{\textrm{e}^{{\left( {\mathfrak {A}-2}\right) }\beta }-1}{\mathfrak {A}-2}, \end{aligned}$$
(32)
$$\begin{aligned} \beta ^*&= \frac{\ln {\left( {\mathfrak {A}-1}\right) }}{\mathfrak {A}-2}, \end{aligned}$$
(33)
$$\begin{aligned} {{\varvec{\varGamma }}^{\prime *}}&= \mathfrak {A}-1. \end{aligned}$$
(34)

The asymptotic behavior for large \(\mathfrak {A}\) is the same as for the complete graph. For \(\mathfrak {A}= 2\), the expressions are ill-defined, but the correct expressions coincide with the limits:

$$\begin{aligned} {\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }&= \ln \beta , \end{aligned}$$
(35)
$$\begin{aligned} \beta ^*&= 1, \end{aligned}$$
(36)
$$\begin{aligned} {{\varvec{\varGamma }}^{\prime *}}&= 1. \end{aligned}$$
(37)

In the biallelic case \(\mathfrak {A}= 2\) this describes accessibility percolation on the directed hypercube, which was considered by Hegarty and Martinsson (2014). In this case \(\beta ^*= 1\), which implies that the directed hypercube is marginally accessible under the HoC model (Franke et al. 2011). For the biallelic case not only the critical value, but also the leading order corrections in the threshold function are known (Hegarty and Martinsson 2014) and coincide with the order \(\frac{\ln L}{L}\) contribution in our candidate threshold function and the value of \({{\varvec{\varGamma }}^{\prime *}}\) given above.

3.3 Path graph

The complete graph is in some sense the best-case scenario for accessibility. On the opposite side of the spectrum of possible (undirected) allele graphs one can choose the path graph on \(\mathfrak {A}\) vertices, shown at the bottom of Fig. 3. In this case the distance between the two end points increases linearly with the number of alleles and there is a unique order in which mutations on a locus must be applied. This causes accessibility to become very low. For \(\mathfrak {A}= 2\) the path graph is identical to the complete graph. However, already for \(\mathfrak {A}= 3\) we find

$$\begin{aligned} \beta ^*= \frac{\ln {\left( {3+2\sqrt{2}}\right) }}{\sqrt{2}} \approx 1.25. \end{aligned}$$
(38)

Since \(\beta ^*\) represents a fitness quantile which must lie between 0 and 1, this value implies that the path graph on three vertices can never be accessible for any fitness difference if (almost) all loci need to mutate from one end of the graph to the other. For higher \(\mathfrak {A}\) this effect becomes more pronounced. As a possible biological application of the path graph the description of copy-number variants of genes can be mentioned (Altenberg 2015).

Since the complete graph on two vertices without reversions has \(\beta ^*= 1\) as shown before in (33) and adding edges can only decrease \(\beta ^*\), it is actually required that the distance between \({a_{l}}\) and \({b_{l}}\) on the allele graph is at least 2 in order for \(\beta ^*> 1\) to be possible.

3.4 Example of non-trivial irregular type

Many graphs seem to be of completely regular type in the sense that no matter which sequence of pairs \(({{\varvec{a}}_L},{{\varvec{b}}_L})\) are chosen, the problem is always of regular type. Martinsson (2018) considered different sufficient conditions on graphs to have this property. But he also lists the smallest graph, of order 4, which does not have it. While this example demonstrates that it is possible to have problems of irregular type, it can also be used to generate semi-regular, but not regular, problem types by carefully interpolating the divergence matrix (2) between a regular and irregular type pair of alleles.

However, the example shown by Martinsson turns out to have \(\beta ^*> 1\), which automatically implies asymptotic inaccessibility in the accessibility percolation context due to the defined range of \(\beta = {F_{{{\varvec{b}}_L}}}-{F_{{{\varvec{a}}_L}}}\) as a difference of uniform random variables. We therefore searched for the smallest graph without the regularity property and \(\beta ^*\le 1\) numerically and found the example in Fig. 5 which has \(\beta ^*\approx 0.983\).

Fig. 5
figure 5

Example of an allele graph that leads to a problem setup of irregular type with \(\beta ^*< 1\)

Fig. 6
figure 6

Allele graph constructed from possible point-mutations on codons. Two amino acids are connected by an arrow iff there is a possible point mutation on a single nucleotide that changes one into the other

3.5 Genetic code

While the complete graph with \(\mathfrak {A}= 4\) may serve as a model for the allele graph of single-nucleotide mutations on DNA or RNA, the expected effect of such a substitution depends significantly on whether or not it changes the amino acid that is encoded by the corresponding three-nucleotide codon. Mutations not affecting the encoded amino acids are known as synonymous. To specifically model the fitness effects of non-synonymous point mutations we therefore consider the allele graph of all amino acids with edges representing the mutual reachability by single-nucleotide substitutions (Fig. 6).

This graph is considerably less symmetric than the complete graph and in particular the resulting quantities \(\beta ^*\) and \({{\varvec{\varGamma }}^{\prime *}}\) will depend on the particular choices of the path endpoints \({{\varvec{a}}_L}\) and \({{\varvec{b}}_L}\) rather than simply on their distance. We consider here all pairings of amino acids \({a_{l}}\) and \({b_{l}}\), assuming them to be equal for all loci. Other cases may be interpolated from these. The results are shown in Table 2. Whether the given values determine the asymptotic behavior of accessibility exactly depends on whether the regularity criteria relating to Martinsson’s function (9) are satisfied. Due to the degree of the graph we limited ourselves to numerical tests, which did not indicate any violation of the criteria, although such violations may be more subtle than our tests could verify.

The critical point \(\beta ^*\) and in particular the expected walk length \(\beta ^*{{\varvec{\varGamma }}^{\prime *}}L\) are, as one would expect, strongly correlated with the distance between alleles. The only distance-3 pair of amino acids is Tyr/Met which also corresponds to the largest walk length with a value of \(\beta ^*{{\varvec{\varGamma }}^{\prime *}}\approx 4.7567\). All other amino acids lie at mutual distance 1 or 2. Nonetheless, the critical point \(\beta ^*\approx 0.4527\) for the distance-3 pair lies slightly below that of the distance-2 pair Asp/Met with \(\beta ^*\approx 0.4570\), demonstrating that the overall structure of the allele graph can have a significant impact on accessibility beyond distance.

For comparison, the accessibility of paths between any pair of codons can be obtained from the values for the complete graph with 4 alleles (Table 1). This gives \(\beta ^*\approx 0.51\), while accounting for the multiplication of three bases per codon yields the expected mean critical walk length per codon as \(\beta ^*{{\varvec{\varGamma }}^{\prime *}}\approx 5.5\). Relative to the complete graph on four alleles the codon graph is allowing arbitrary synonymous mutations without cost, causing the reduction in the critical fitness difference as well as the length of walks at the critical point. However, when compared to the complete graph on 21 alleles with \(\beta ^*\approx 0.154\) and critical walk length factor 3.22 (Table 1), which would permit direct mutations between arbitrary amino acids, the fitness cost is still significantly higher on the codon graph, whereas the critical walk length on the codon graph is scattered around that of the 21-allele complete graph.

Table 2 Top: critical fitness difference \(\beta ^*\) for accessible paths between homopolymer amino acid sequences consisting of the indicated pairs

4 First moment bound and walk length

We start with a proof for the first part of Theorem 1 from an upper bound for accessibility based on the mean number of accessible paths, or rather the mean number of quasi-accessible walks. We define the term quasi-accessible as a generalization of the notion of accessibility used up to now as explained in the following.

4.1 Quasi-accessibility

In the original definition of accessibility, a non-self-avoiding walk is never accessible, because it would have to visit the same fitness value twice, which makes it impossible for the walk to have strictly increasing fitness. Handling self-avoidance is non-trivial. To remedy this in a simpler manner, instead of considering self-avoiding paths on \({\mathcal {A}^L}\), we consider an extension of \({\mathcal {A}^L}\) to \({\mathcal {A}^L}'\) as follows:

Definition 14

The extended genotype space of a genotype space \({\mathcal {A}^L}\) is the simple directed graph \({\mathcal {A}^L}'\) with vertex set \({\mathbb {A}^L}' = {\mathbb {A}^L}\times \mathbb {N}\) and an arrow from \(({\varvec{v}},n)\) to \(({\varvec{w}},m)\) iff there is an arrow from \({\varvec{v}}\) to \({\varvec{w}}\) in \({\mathcal {A}^L}\).

In other words we duplicate every genotype a countable infinite number of times in such a way that traversal of one of its copies can always be replaced by traversal of another copy. The 1-section containing all vertices of the form \(({\varvec{v}}, 1)\) can be identified with the vertices on \({\mathcal {A}^L}\). We then assign each of the vertices in \({\mathcal {A}^L}'\) i.i.d. fitness values:

Definition 15

An extended HoC model over a genotype space \({\mathcal {A}^L}\) is the HoC model over the extended genotype space of the genotype space \({\mathcal {A}^L}\).

The mentioned 1-section then corresponds to the original HoC model and we can identify realizations of the extended HoC model with the corresponding realizations of the original HoC model with equal fitness values on the 1-section. All other fitness values do not affect this underlying model. However, it is convenient to introduce these additional fitness values for the following reasons.

We define the following map of walks on \({\mathcal {A}^L}\) to \({\mathcal {A}^L}'\). Each self-avoiding walk is mapped to the corresponding walk on the 1-section of \({\mathcal {A}^L}'\). But instead of mapping non-selfavoiding walks from \({\mathcal {A}^L}\) to the 1-section of \({\mathcal {A}^L}'\), we can make use of the additional vertex copies to replace all vertices that are visited multiple times in \({\mathcal {A}^L}\) with distinct copies in \({\mathcal {A}^L}'\). To make this unique, we assume that the n-th visit of vertex \({\varvec{v}}\) in \({\mathcal {A}^L}\) is mapped to the vertex \(({\varvec{v}},n)\) in \({\mathcal {A}^L}'\), except if \({\varvec{v}}\) is the final vertex of the walk, in which case we map the n-th visit in reverse order to \(({v_{}},n)\) in \({\mathcal {A}^L}'\). The resulting walk is always selfavoiding in \({\mathcal {A}^L}'\) and the special case assures that every walk in \({\mathcal {A}^L}\) is mapped to a walk with endpoints on the 1-section in \({\mathcal {A}^L}'\).

Definition 16

A walk on \({\mathcal {A}^L}\) is quasi-accessible if the corresponding mapped walk on \({\mathcal {A}^L}'\) per the rules above is accessible on \({\mathcal {A}^L}'\). A genotype is said to be quasi-accessible from another if there exists a quasi-accessible walk from the latter to the former.

Definition 17

A walk on \({\mathcal {A}^L}'\) is \(valid \) if there is a walk in \({\mathcal {A}^L}\) which is mapped to it according to the rules above. An (extended) genotype \({\varvec{h}}'\) on \({\mathcal {A}^L}'\) is said to be valid-accessible from a genotype \({\varvec{g}}'\) on \({\mathcal {A}^L}'\) if there exists a valid-accessible walk from the latter to the former on \({\mathcal {A}^L}'\).

With these definitions the probability of any walk on \({\mathcal {A}^L}\), whether self-avoiding or not, to be quasi-accessible is the same, only depending on the length of the walk. Furthermore the notion of valid-accessibility on \({\mathcal {A}^L}'\) coincides with quasi-accessibility on \({\mathcal {A}^L}\). We denote the number of valid-accessible walks from \(({\varvec{v}},1) \in {\mathbb {A}^L}'\) to \(({\varvec{w}},1) \in {\mathbb {A}^L}'\), or equivalently quasi-accessible walks from \({\varvec{v}}\) to \({\varvec{w}}\) with \({\varvec{v}},{\varvec{w}}\in {\mathbb {A}^L}\), by \(\tilde{{\varvec{Z}}}_{{{\varvec{v}}}{{\varvec{w}}}}\). Additionally, while quasi-accessibility is different from accessibility for individual walks on \({\mathcal {A}^L}\), accessibility and quasi-accessibility of one genotype from another on \({\mathcal {A}^L}\) coincide:

Lemma 1

$$\begin{aligned} {\mathbb {P}\left[ {\varvec{Z}}_{{{\varvec{v}}}{{\varvec{w}}}} \ge 1\right] } = {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{v}}}{{\varvec{w}}}} \ge 1\right] }. \end{aligned}$$
(39)

Proof

If \({\varvec{w}}\) is accessible from \({\varvec{v}}\), then there exists a walk from \({\varvec{v}}\) to \({\varvec{w}}\) which is accessible and therefore also quasi-accessible, implying

$$\begin{aligned} {\mathbb {P}\left[ {\varvec{Z}}_{{{\varvec{v}}}{{\varvec{w}}}} \ge 1\right] } \le {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{v}}}{{\varvec{w}}}} \ge 1\right] }. \end{aligned}$$
(40)

If \({\varvec{w}}\) is quasi-accessible from \({\varvec{v}}\), then there exists a valid walk on \({\mathcal {A}^L}'\) from \(({\varvec{v}},1)\) to \(({\varvec{w}},1)\) which is accessible (with respect to \({\mathcal {A}^L}'\)). Removing all vertices \(({\varvec{u}},n)\) with \(n \ne 1\) from this walk results in another walk on \({\mathcal {A}^L}'\) which is completely located on the 1-section due to the validity requirement. Because this walk is obtained by only removal of vertices, it is also valid-accessible and because it must be self-avoiding on the 1-section, it is also an accessible walk from \({\varvec{v}}\) to \({\varvec{w}}\) implying the other side of the equality

$$\begin{aligned} {\mathbb {P}\left[ {\varvec{Z}}_{{{\varvec{v}}}{{\varvec{w}}}} \ge 1\right] } \ge {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{v}}}{{\varvec{w}}}} \ge 1\right] }. \end{aligned}$$
(41)

\(\square \)

This implies that we can restrict our investigation to quasi-accessibility.

4.2 Proof of the lower side condition of Theorem 1.

In order to give an upper bound on \(\beta \)-quasi-accessibility and with that a proof of the first part of Theorem 1, we will consider the mean number of quasi-accessible walks from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\). Each walk of length \(N\) from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\) on \({\mathcal {A}^L}\) is \(\beta \)-quasi-accessible with probability

$$\begin{aligned} \frac{\beta ^{N-1}}{(N-1)!} \end{aligned}$$

where the numerator accounts for the probability that all inner vertices of the walk are found inside the range of fitness values \({F_{{{\varvec{a}}_L}}}\) to \({F_{{{\varvec{b}}_L}}}\) and the denominator accounts for the increasing order required on these values. The number of walks taking \({n_{}}\) steps from \({a_{l}}\) to \({b_{l}}\) on one locus \(l\) is given by \({\left( A^{{n_{}}}\right) _{{{a_{l}}}{{b_{l}}}}}\). A walk of length \(N\) could take each step on any of the loci, so that the total number of walks of length \(N\) can be written as

$$\begin{aligned} \sum _{{n_{1}}+\cdots +{n_{L}} = N} \left( {\begin{array}{c}N\\ {n_{1}},\ldots ,{n_{L}}\end{array}}\right) \prod _{l=1}^L{\left( A^{{n_{l}}}\right) _{{{a_{l}}}{{b_{l}}}}} \end{aligned}$$
(42)

where \(\left( {\begin{array}{c}N\\ {n_{1}},\ldots ,{n_{L}}\end{array}}\right) \) is the multinomial coefficient accounting for the different orderings of steps on individual loci. Multiplication of this expression with the probability of quasi-accessibility of each such walk gives the mean number of quasi-accessible paths

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } = \sum _{N=0}^\infty \sum _{{n_{1}}+\cdots +{n_{L}} = N} \left( {\begin{array}{c}N\\ {n_{1}},\ldots ,{n_{L}}\end{array}}\right) \frac{\beta ^{N-1}}{{\left( {N-1}\right) }!} \prod _{l=1}^L{\left( A^{{n_{l}}}\right) _{{{a_{l}}}{{b_{l}}}}}. \end{aligned}$$
(43)

The term \({\left( {N-1}\right) }!\) can be reduced to \(N!\) by introduction of a derivative

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } = \partial _\beta \sum _{N=0}^\infty \sum _{{n_{1}}+\cdots +{n_{L}} = N} \left( {\begin{array}{c}N\\ {n_{1}},\ldots ,{n_{L}}\end{array}}\right) \frac{\beta ^{N}}{N!} \prod _{l=1}^L{\left( A^{{n_{l}}}\right) _{{{a_{l}}}{{b_{l}}}}} \end{aligned}$$
(44)

and redistributing all the factorials and \(\beta ^N\) into the product yields

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } = \partial _\beta \sum _{N=0}^\infty \sum _{{n_{1}}+\cdots +{n_{L}} = N} \prod _{l=1}^L\frac{\beta ^{{n_{l}}}}{{n_{l}}!}{\left( A^{{n_{l}}}\right) _{{{a_{l}}}{{b_{l}}}}}. \end{aligned}$$
(45)

Finally the sums and the product can be interchanged and

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }&= \partial _\beta \prod _{l=1}^L\sum _{{n_{}}=0}^\infty \frac{\beta ^{{n_{}}}}{{n_{}}!}{\left( A^{{n_{}}}\right) _{{{a_{l}}}{{b_{l}}}}} = \partial _{\beta }\prod _{l=1}^L{\left( \textrm{e}^{{\beta }A}\right) _{{{a_{l}}}{{b_{l}}}}} = \partial _\beta \textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }} \end{aligned}$$
(46)
$$\begin{aligned}&= L{\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) } \textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }}. \end{aligned}$$
(47)

We now choose \(\beta \) based on Theorem 1 as

$$\begin{aligned} \beta = {\hat{\beta }}- {{\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) }}^{-1}\frac{\ln L}{L} - \frac{|g(L)|}{L} \end{aligned}$$
(48)

where \(g(L) = \omega (1)\). Because \({\hat{\beta }}\) converges to a positive value per Proposition 4, around which \(\varGamma _{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\) and its derivatives are bounded, a Taylor expansion of \({\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\) gives

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } = L{\left( {{\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) } + \mathcal {O}{{\left( {\beta -{\hat{\beta }}}\right) }}}\right) } \textrm{e}^{L{\left( {{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) } + {\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) }{\left( {\beta -{\hat{\beta }}}\right) } + \mathcal {O}{\left( {{{\left( {\beta -{\hat{\beta }}}\right) }^2}}\right) } }\right) } }. \end{aligned}$$
(49)

By definition \({\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) } = 0\) and by Proposition 4\({\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) }\) converges to a positive value \({{\varvec{\varGamma }}^{\prime *}}\), so that inserting the difference \(\beta -{\hat{\beta }}\) then gives

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } = ({{\varvec{\varGamma }}^{\prime *}}+o(1))\textrm{e}^{-{{\varvec{\varGamma }}^{\prime *}}|g(L)|}, \end{aligned}$$
(50)

implying with \(g(L) = \omega (1)\) that

$$\begin{aligned} \lim _{L\rightarrow \infty }{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } = 0. \end{aligned}$$
(51)

Then, by Markov’s inequality we have \({\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1\right] } \le {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }\), implying that quasi-accessibility converges to zero, which then by Lemma 1 also implies that accessibility of \({{\varvec{b}}_L}\) from \({{\varvec{a}}_L}\) converges to zero, proving part 1 of Theorem 1. \(\square \)

Corollary 1

By replacing \(-|g(L)|\) by a constant value \(\eta \) in the above, the expected number of quasi-accessible walks will converge to a non-zero value and increasing \(\eta \) allows to arbitrarily increase the limit.

4.3 Proof of Theorem 3

A more refined version of the previous argument can be used to prove Theorem 3. Specifically the expected number of \(\beta \)-quasi-accessible walks can be separated into intervals of walk lengths. Let \(h_N\) be the expected number of \(\beta \)-quasi-accessible walks of length N with

$$\begin{aligned} \beta = {\hat{\beta }}- {{\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{{\hat{\beta }}}}\right) }}^{-1}\frac{\ln L}{L} + \frac{\eta }{L} \end{aligned}$$
(52)

as in Theorem 3 for some fixed value \(\eta \). This number is an expectation value over realizations of fitness values, but in the following we will consider it as just a number indexed by some number \(N\) representing the walk length. Summation of all of these numbers then yields the total expected number of \(\beta \)-quasi-accessible walks which we calculated already above:

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } = \partial _{\beta } \textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }} = \sum _{N=1}^\infty h_N. \end{aligned}$$
(53)

We can interpret this as the value \(\phi (1)\) of the function

$$\begin{aligned} \phi (z) = \partial _{z\beta }\textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{z\beta }}\right) }} = \sum _{N=1}^\infty h_Nz^{N-1}. \end{aligned}$$
(54)

This function can be viewed as an (ordinary) generating function for the sequence \(h_N\) shifted by one. The generating function here is not related to the probability distribution of fitness values, but is rather to be understood as simply a counting tool that separates the total expectation value into slots for different walks lengths using the additivity of the expectation value.

The effect of the derivative \(\partial _{z\beta }\) in the generating function can be reversed by integration of each of the monomials, so that

$$\begin{aligned} {\tilde{\phi }}(z) = \textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{z\beta }}\right) }} = \sum _{N=1}^\infty \frac{\beta }{N} h_Nz^{N} = \sum _{N=1}^\infty {\tilde{h}}_Nz^N\end{aligned}$$
(55)

is the generating function of the unshifted \(h_N\) multiplied by \(\frac{\beta }{N}\), which is another sequence that we define as \({\tilde{h}}_N\). Normalizing \({\tilde{\phi }}(z)\) through division by \({\tilde{\phi }}(1) = \textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }}\) turns the generating function into a probability generating function over the parameter \(N\) as random variable and this allows us to apply theorems from probability theory. Again, this probability is not related to the distribution of fitness values, but is introduced here artificially as a counting tool. The integrated (probability) generating function factorizes over loci as

$$\begin{aligned} \frac{{\tilde{\phi }}(z)}{{\tilde{\phi }}(1)} = \prod _{l=1}^L\frac{ \textrm{e}^{\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{z\beta }}\right) }} }{ \textrm{e}^{\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{\beta }}\right) }} }. \end{aligned}$$
(56)

As a consequence the random variable \(N\) under the generating function’s distribution can be written as a sum \(N= \sum _{l=1}^L{n_{l}}\), where \({n_{l}}\) are independent random variables with probability generating functions

$$\begin{aligned} \textrm{e}^{\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{z\beta }}\right) } -\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{\beta }}\right) }} = \frac{ \sum _{{n_{l}}=0}^\infty \frac{\beta ^{{n_{l}}}}{{n_{l}}!}{\left( A^{{n_{l}}}\right) _{{{a_{l}}}{{b_{l}}}}} z^{{n_{l}}} }{ \sum _{{n_{l}}=0}^\infty \frac{\beta ^{{n_{l}}}}{{n_{l}}!}{\left( A^{{n_{l}}}\right) _{{{a_{l}}}{{b_{l}}}}} }. \end{aligned}$$
(57)

This can be seen going in the reverse direction as the generating function of the sum of independent random variables is the product of the individual generating functions of the summands. Because the degree of \(A\) is bounded by \(\varDelta \), \({\left( A^{{n_{l}}}\right) _{{{a_{l}}}{{b_{l}}}}} \le \varDelta ^{{n_{l}}}\) and the tail of the distribution is dominated by an exponential. This bound is also independent of the chosen loci \({a_{l}}\) and \({b_{l}}\) and with the chosen \(\beta \) converging in \(L\), the central limit theorem applies to the sum \(N\). The mean \(L\mu \) and variance \(L\sigma ^2\) of \(N\) under this distribution can be obtained from the first and second derivatives of the probability generating function as

$$\begin{aligned} \mu&= \beta {\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) } \end{aligned}$$
(58)
$$\begin{aligned} \sigma ^2&= \beta {\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) } + {\beta }^2{\varvec{\varGamma }}''_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }, \end{aligned}$$
(59)

and the central limit theorem implies that for constants \(c>0\):

$$\begin{aligned} \sum _{|N-\mu L|\ge c\sigma \sqrt{L}} {\tilde{h}}_N= 2\textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }}\varPhi (-c){\left( {1+o(1)}\right) } \end{aligned}$$
(60)

where \(\varPhi \) is the standard normal CDF. Since the sum’s upper bound is asymptotic to \(\mu L\), for all \({\tilde{h}}_N\) terms appearing in the sum \({\tilde{h}}_N\ge \frac{\beta }{\mu L}h_N\), so that

$$\begin{aligned} \sum _{|N-\mu L|\ge c\sigma \sqrt{L}} h_N\le 2\frac{\mu L}{\beta }\textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }}\varPhi (-c){\left( {1+o(1)}\right) }. \end{aligned}$$
(61)

In particular with the choice of the threshold function of Theorem 3 for \(\beta \):

$$\begin{aligned} \sum _{|N-\mu L|\ge c\sigma \sqrt{L}} h_N \textrm{e}^{\eta }\le 2\eta \varPhi (-c){\left( {1+o(1)}\right) }. \end{aligned}$$
(62)

This allows one to reduce the mean number of \(\beta \)-quasi-accessible walks of length outside the interval \(\mu L\pm c\sigma \sqrt{L}\) to any arbitrarily small value by choosing c large enough. In other words, if there are \(\beta \)-quasi-accessible walks at the suggested threshold function, then they are of length \(\mu L\) with fluctuations of at most order of \(\sqrt{L}\). Since all \(\beta \)-accessible walks are also \(\beta \)-quasi-accessible walks, the same length constraint then also applies to accessible paths at the threshold function, completing the proof of Theorem 3. \(\Box \)

5 Proof of Theorem 2

The upper bound on accessibility obtained from the expected value does not take into account any dependence between walks. We can improve the bound by including some of the dependencies. This will make it possible to prove our Theorem 2, i.e. show that \(\beta _\text {c}> \beta ^*\) for irregular types.

Let in this section \(0< r< 1\) and \(0< s< 1\) be constant. The intention is to choose them later such that \({\mathfrak {M}^*({s},{r},{\beta ^*})} > 0\) as application to the irregular type. This choice is always possible in the irregular case since by continuity \({\mathfrak {M}^*({s},{r},{\beta ^*})}\) cannot be strictly positive only on the boundaries.

Recall that the parameters \(r\) and \(s\) determine the fitness spanned by three segments of each walk with \(s\) determining the fitness fraction spanned by the middle segment, and \(r\) determining the distribution of the remaining fitness span onto the first and last segment. More concretely the intended fitness span of the first segment is \(\beta {{\bar{s}}}r\), of the second \(\beta s\) and the third \(\beta {{\bar{s}}}{{\bar{r}}}\), adding up to the full fitness span \(\beta \) that needs to be crossed [see Sect. 2.2, Eq. (14)].

For each arrow on the genotype space we can consider the interval formed by the fitness values of the two genotypes incident to it. If a walk from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\) is accessible, then it contains exactly one arrow with a fitness interval containing the fitness value \(\beta {{\bar{s}}}r\). Let this arrow be \(({\varvec{x}},{\varvec{x}}')\) and let \(S^{12}_{{\varvec{x}},{\varvec{x}}'}\) be an indicator variable for this fitness value falling on to the arrow \(({\varvec{x}}, {\varvec{x}}')\). Similarly there is exactly one arrow containing the fitness value \(\beta {\left( {1-{{\bar{s}}}{{\bar{r}}}}\right) }\). Let this arrow be \(({\varvec{y}},{\varvec{y}}')\) and the corresponding indicator \(S^{23}_{{\varvec{y}},{\varvec{y}}'}\). These two arrows segment the walk in the closest possible way according to the intended fitness spans mentioned above. A walk is accessible only if each of the three segments \({{\varvec{a}}_L}\rightarrow {\varvec{x}}\), \({\varvec{x}}' \rightarrow {\varvec{y}}\) and \({\varvec{y}}' \rightarrow {{\varvec{b}}_L}\) are accessible. In the following we refer to these segments as segment 1, 2 and 3 respectively. To obtain an upper bound on the accessibility of \({{\varvec{b}}_L}\) from \({{\varvec{a}}_L}\) it is therefore sufficient to form a union bound of the form:

$$\begin{aligned}&{\mathbb {P}\left[ {\varvec{Z}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1\right] } \end{aligned}$$
(63)
$$\begin{aligned} \le&\sum _{({\varvec{x}},{\varvec{x}}'),({\varvec{y}},{\varvec{y}}')\in {\mathcal {A}^L}} {\mathbb {P}\left[ S^{12}_{{\varvec{x}},{\varvec{x}}'} \wedge S^{23}_{{\varvec{y}},{\varvec{y}}'} \wedge {\varvec{Z}}_{{{{\varvec{a}}_L}}{{\varvec{x}}}} \ge 1 \wedge {\varvec{Z}}_{{{\varvec{x}}'}{{\varvec{y}}}} \ge 1 \wedge {\varvec{Z}}_{{{\varvec{y}}'}{{{\varvec{b}}_L}}} \ge 1\right] }. \end{aligned}$$
(64)

We can now separate the expectation over the fitness of the intermediate genotypes \({\varvec{x}}\), \({\varvec{x}}'\), \({\varvec{y}}\) and \({\varvec{y}}'\),

$$\begin{aligned}&{\mathbb {P}\left[ {\varvec{Z}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1\right] } \end{aligned}$$
(65)
$$\begin{aligned}&\quad \le \sum _{({\varvec{x}},{\varvec{x}}'),({\varvec{y}},{\varvec{y}}')\in {\mathcal {A}^L}} {\mathbb {E}\left[ {\mathbb {P}\left[ {\varvec{Z}}_{{{{\varvec{a}}_L}}{{\varvec{x}}}} \ge 1 \wedge {\varvec{Z}}_{{{\varvec{x}}'}{{\varvec{y}}}} \ge 1 \wedge {\varvec{Z}}_{{{\varvec{y}}'}{{{\varvec{b}}_L}}} \ge 1 | \cdot \right] }\right] }, \end{aligned}$$
(66)

where the dot in the probability indicates conditioning on the fitness values at the end points. The expectation is over all such fitness values satisfying the conditions \(S^{12}_{{\varvec{x}},{\varvec{x}}'}\) and \(S^{23}_{{\varvec{y}},{\varvec{y}}'}\).

As a result of this conditioning, the quasi-accessibilities of the three segments mentioned in the equation are negatively dependent, so that the joint probability of the events can be upper bounded by the product of individual probabilities. In order for a segment to be accessible under the conditioning of the fitness values at the end points, all internal fitness values on the segmnent must fall into the fitness range between the end points and the internal fitness values must be increasingly ordered. Because the fitness values are i.i.d. in the HoC model these two properties are independent. Now suppose we condition on one (or two) of the segments being accessible with any particular choice of accessible walks. This is equivalent to conditioning all the internal fitness values of these walks in an appropriate manner as well. The accessibility requires these internal fitness values to be constrained to the fitness range of the segment’s end points, which makes them unavailable as internal vertices of the remaining segment(s). However, all other fitness values are i.i.d. and unaffected by this conditioning. The probability that the remaining segment(s) is (are) then accessible is therefore smaller than if no conditioning of the internal vertices of the other accessible segments had been applied, effectively only removing walks through accessible vertices of the conditioned segments from the set of candidate walks. For example if \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{\varvec{x}}}}\) is at least 1, then \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}'}{{\varvec{y}}}} \ge 1\) becomes less likely since the existence of a quasi-accessible walk from \({{\varvec{a}}_L}\) to \({\varvec{x}}\) implies that some fitness values of other genotypes fall in the range \([{F_{{{\varvec{a}}_L}}},{F_{{\varvec{x}}}}]\), excluding them for consideration in the range \([{F_{{\varvec{x}}'}}, {F_{{\varvec{y}}}}]\) required for them to be part of a quasi-accessible walk from \({\varvec{x}}'\) to \({\varvec{y}}\). Consequently:

$$\begin{aligned} {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1\right] } \le \sum _{{\left( {{\varvec{x}},{\varvec{x}}'}\right) },{\left( {{\varvec{y}},{\varvec{y}}'}\right) }} {\mathbb {E}\left[ {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{\varvec{x}}}} \ge 1 | \cdot \right] } {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}'}{{\varvec{y}}}} \ge 1 | \cdot \right] } {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{y}}'}{{{\varvec{b}}_L}}} \ge 1 | \cdot \right] } \right] }. \end{aligned}$$
(67)

Since probabilities lie in [0, 1], we can use the upper bound \(x \le x^{1-\alpha }\) for \(0< \alpha < 1\) on the middle factor and afterwards we can apply Markov’s inequality to all three terms to obtain

$$\begin{aligned} {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1\right] } \le \sum _{{\left( {{\varvec{x}},{\varvec{x}}'}\right) },{\left( {{\varvec{y}},{\varvec{y}}'}\right) }} {\mathbb {E}\left[ {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{\varvec{x}}}} \ge 1 | \cdot \right] } {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}'}{{\varvec{y}}}} \ge 1 | \cdot \right] }^{1-\alpha } {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{y}}'}{{{\varvec{b}}_L}}} \ge 1 | \cdot \right] } \right] }. \end{aligned}$$
(68)

The remaining inner expectation values depend only on the differences of the fitness values that they are conditioned on, not the actual placement of that difference. We introduce the following quantities:

$$\begin{aligned} \epsilon _{1}&= \beta {{\bar{s}}}r- {\left( {{F_{{\varvec{x}}}}-{F_{{{\varvec{a}}_L}}}}\right) }, \end{aligned}$$
(69)
$$\begin{aligned} \epsilon _{2}&= \beta s- {\left( {{F_{{\varvec{y}}}}-{F_{{\varvec{x}}'}}}\right) }, \end{aligned}$$
(70)
$$\begin{aligned} \epsilon _{3}&= \beta {{\bar{s}}}{{\bar{r}}}- {\left( {{F_{{{\varvec{b}}_L}}}-{F_{{\varvec{y}}'}}}\right) }. \end{aligned}$$
(71)

These quantities measure how much the fitness difference allocated to one of the three walk segments differs from what it would be assigned if \(r\) and \(s\) determined it exactly. For example the \(({\varvec{x}},{\varvec{x}}')\) edge is required to contain the fitness value \(\beta {{\bar{s}}}r\). Therefore the first walk segment can span a fitness distance of at most \(\beta {{\bar{s}}}r\), but this happens exactly only if \({F_{{\varvec{x}}}} - {F_{{{\varvec{a}}_L}}} = \beta {{\bar{s}}}r\) is chosen. All other valid choices set the fitness value lower than this and \(\epsilon _1\) measures the reduction of the segment’s length. As it will turn out only the point with all \(\epsilon _i\) equal to zero contributes to the expectation value in leading order. Intuitively any constant offset from the intended segment length corresponds to an effective reduction of \(\beta \) by a constant, resulting in an exponentially lower likelihood of walks being quasi-accessible. Nonetheless we will carry the \(\epsilon _i\) through the calculation.

The remaining inner expectation values are of the same form as the simple expectation of walks from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\) calculated in the previous section:

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{v}}}{{\varvec{w}}}}|t\right] } = \partial _t\textrm{e}^{L{\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }} = {\varvec{\varGamma }}'_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }L\textrm{e}^{L{\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }}. \end{aligned}$$
(72)

The expectation values are dominated by the exponential terms \(\textrm{e}^{L{\varvec{\varGamma }}_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }}\) with an additional linear factor \(L\) resulting from the derivative. However, the derivative also adds the term \({\varvec{\varGamma }}'_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) }\). As an average over loci it can be seen that pointwise in \(t\) and uniformly over \({\varvec{v}}\) and \({\varvec{w}}\), this quantity is bounded by a constant from above. However the bound is not uniform in \(t\). At \(t= 0\) it diverges, as can be seen from the expansion

$$\begin{aligned} L{\varvec{\varGamma }}'_{{{\varvec{v}}}{{\varvec{w}}}}{\left( {{t}}\right) } \sim \frac{{{\varvec{d}}_{{{\varvec{v}}}{{\varvec{w}}}}}}{t}. \end{aligned}$$
(73)

To avoid this issue we rewrite the expectation value including the sum resulting from application of the product rule of differentiation:

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{v}}}{{\varvec{w}}}}|t\right] } = \sum _{l'=1}^L{\left( A\textrm{e}^{{t}A}\right) _{{{v_{l'}}}{{w_{l'}}}}} \prod _{l\ne l'}\textrm{e}^{\varGamma _{{{v_{l}}}{{w_{l}}}}{\left( {{t}}\right) }}. \end{aligned}$$
(74)

In each summand the value is a product over terms, each of which depends only on quantities on a single locus and the bulk of the contributions of loci contribute simply the exponential \(\textrm{e}^{\varGamma _{{{v_{l}}}{{w_{l}}}}{\left( {{t}}\right) }} = {\left( \textrm{e}^{{t}A}\right) _{{{v_{l}}}{{w_{l}}}}}\). Only the locus \(l'\) gives a different contribution, namely the derivative of the exponential term, \({\left( A\textrm{e}^{{t}A}\right) _{{{v_{l'}}}{{w_{l'}}}}}\).

Our goal is to bring Eq. (68) into the form of a sum over products, such that the product factorizes in the same sense as it does for a single expectation value. In particular the current form is a sum of a product of three expectation values. If we expand each expectation as shown in Eq. (74), we obtain three sums, each accounting for one special locus on which the corresponding derivative is taken. We name these special loci \(l_1\), \(l_2\) and \(l_3\), corresponding to the means in Eq. (68) in the order they appear there. The sum in the middle term can be taken out of the \({\left( {\cdot }\right) }^{1-\alpha }\) form to give an upper bound, because \(1-\alpha \in [0,1]\) and therefore the form is subadditive. Having done so, the sum over the pair of edges on the genotype space may similarly be factorized over loci. Each edge on the genotype graph corresponds to a step on one locus. Therefore it is sufficient to sum over individual genotypes together with another special locus, and one edge on the allele graph corresponding to that locus. We denote the sum over loci for these two edges \(l_{12}\) and \(l_{23}\) respectively. The initial sum then factorizes over loci:

$$\begin{aligned} {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1\right] } \le \sum _{l_1,l_2,l_3,l_{12},l_{23} = 1}^L{\mathbb {E}\left[ \prod _{l=1}^L\mathfrak {F}_{l} \right] }. \end{aligned}$$
(75)

Here \(\mathfrak {F}_{l}\) is the resulting factor collecting all sums over quantities on locus \(l\) and all factors of the product of the three expectations that are functions of quantities on locus \(l\), as well as potentially e.g. a form \(\sum _{{x_{l}}'}{\left( A^{}\right) _{{{x_{l}}}{{x_{l}}'}}}\) if \(l= l_{12}\). \(\mathfrak {F}_{l}\) is implicitly dependent on \(l_1\), \(l_2\), \(l_3\), since these three variables decide whether the contribution resulting from any of the three expectation values has the usual exponential form or that of its derivative. If \(l\not \in \{l_1,l_2,l_3,l_{12},l_{23}\}\), then the contribution of all three expectation values and the edge sum is of the usual form, i.e. the exponential term of the expectation value and no sum over edges and we give it the name \(\mathfrak {G}_l\):

$$\begin{aligned} \mathfrak {G}_l&= \sum _{{x_{l}},{y_{l}}} \textrm{e}^{ \varGamma _{{{a_{l}}}{{x_{l}}}}{\left( {{\beta {{\bar{s}}}r-\epsilon _{1}}}\right) } +{\left( {1-\alpha }\right) }\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s-\epsilon _{2}}}\right) } +\varGamma _{{{y_{l}}}{{b_{l}}}}{\left( {{\beta {{\bar{s}}}{{\bar{r}}}-\epsilon _{3}}}\right) } } \end{aligned}$$
(76)
$$\begin{aligned}&= \sum _{{x_{l}},{y_{l}}} {\left( \textrm{e}^{{{\left( {\beta {{\bar{s}}}r-\epsilon _{1}}\right) }}A}\right) _{{{a_{l}}}{{x_{l}}}}} {\left( {{\left( \textrm{e}^{{{\left( {\beta s-\epsilon _{2}}\right) }}A}\right) _{{{x_{l}}}{{y_{l}}}}}}\right) }^{1-\alpha } {\left( \textrm{e}^{{{\left( {\beta {{\bar{s}}}{{\bar{r}}}-\epsilon _{3}}\right) }}A}\right) _{{{y_{l}}}{{b_{l}}}}}. \end{aligned}$$
(77)

If \(l\) is equal to any of the set of special loci, then some of these exponential terms will be modified and there might be additional sums. For example if \(l\) is equal to \(l_3\) and \(l_{12}\), but not equal to any of the other special loci, then

$$\begin{aligned} \mathfrak {F}_l= \sum _{{x_{l}},{y_{l}},{x_{l}}'} {\left( \textrm{e}^{{{\left( {\beta {{\bar{s}}}r-\epsilon _{1}}\right) }}A}\right) _{{{a_{l}}}{{x_{l}}}}} {\left( A^{}\right) _{{{x_{l}}}{{x_{l}}'}}} {\left( {{\left( \textrm{e}^{{{\left( {\beta s-\epsilon _{2}}\right) }}A}\right) _{{{x_{l}}'}{{y_{l}}}}}}\right) }^{1-\alpha } {\left( A\textrm{e}^{{{\left( {\beta {{\bar{s}}}{{\bar{r}}}-\epsilon _{3}}\right) }}A}\right) _{{{y_{l}}}{{b_{l}}}}}. \end{aligned}$$
(78)

By assumption \(s\ne 0\) and also \(0 \ne r\ne 1\). Then, with the fixed \(s\) and \(r\), \(\mathfrak {G}_l\) is bounded away from zero everywhere except at the boundary, because as long as one of the matrix exponentials has non-zero argument, it contributes a finite term to the sum by adequate choice of \({x_{l}}\) and \({y_{l}}\) so that the indices of the matrix exponential become \(({a_{l}},{b_{l}})\). More generally \(\mathfrak {G}_l\) is also uniformly bounded over \(s\) and \(r\), since by definition of \(s\) and \(r\) at least one of the matrix exponential arguments must be at least \(\frac{\beta }{3}\), epsilon shifts notwithstanding. This then allows us to write each \(\mathfrak {F}_l\) as a product \(\mathfrak {G}_l \mathfrak {H}_l\), with \(\mathfrak {H}_l\) bounded away from infinity except at the mentioned boundary. In the next section we will use the same approach with a more detailed handling of \(\mathfrak {H}_l\), but here it is sufficient to apply such a simple uniform bound with a constant.

However first we consider the behavior at the boundary where all epsilon shifts force the matrix exponential arguments to become zero. Due to the bounded degree of the graph, as \(t\rightarrow 0\), the diagonal terms of the matrix exponential with argument \(t\) drop to 1 and the off-diagonal ones to 0 uniformly. If \({a_{l}}\ne {b_{l}}\), \(\mathfrak {F}_l\) therefore falls to zero as all the \(\epsilon _i\) reach their maximum boundary and similarly it falls to 1 for \({a_{l}}= {b_{l}}\). Since there is by assumption at least a finite fraction of loci with \({a_{l}}\ne {b_{l}}\), this then implies that eventually, at a finite distance to the boundary

$$\begin{aligned} \prod _{l=1}^L\mathfrak {F}_l \le \mathcal {O}{\left( {{C^L}}\right) } \end{aligned}$$
(79)

for some \(C < 1\). The special loci on which \(\mathfrak {F}_l \ne \mathfrak {G}_l\) are not relevant to this, since there are only finitely many of them and each one is bounded. The contribution to the probability from the boundary is therefore asymptotically zero, since the exponential decay in the integrand cannot be compensated by the additional \(L^5\) factor from the special loci sum.

Returning to the general case away from the boundary, we can bound all \(\mathfrak {H}_l\) with \(l\in \{l_1,l_2,l_3,l_{12},l_{13}\}\) by some constant C uniformly, yielding a factor of at most \(C^5\), while all other \(\mathfrak {H}_l\) are 1. This removes the dependence of the product on the particular choice of the special loci:

$$\begin{aligned} {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1 \right] } \le L^5 C^5 {\mathbb {E}\left[ \prod _{l=1}^L\mathfrak {G}_l\right] } + o(1) = L^5 C^5 {\mathbb {E}\left[ \textrm{e}^{L\mathfrak {T}}\right] } \end{aligned}$$
(80)

with

$$\begin{aligned} \mathfrak {T} = {\left\langle {\ln \sum _{{x_{l}},{y_{l}}}\textrm{e}^{ \varGamma _{{{a_{l}}}{{x_{l}}}}{\left( {{\beta {{\bar{s}}}r-\epsilon _{1}}}\right) } +{\left( {1-\alpha }\right) }\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s-\epsilon _{2}}}\right) } +\varGamma _{{{y_{l}}}{{b_{l}}}}{\left( {{\beta {{\bar{s}}}{{\bar{r}}}-\epsilon _{3}}}\right) } }}\right\rangle }_{l}. \end{aligned}$$
(81)

The \(\epsilon _i\) are always non-negative in the valid domain and \(\mathfrak {T}\) is decreasing in all of them. Therefore we can give an upper bound by setting all of them to 0 and obtain the upper bound on accessibility:

$$\begin{aligned} {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1 \right] } \le \textrm{e}^{L{\left( {\mathfrak {T}_0 + o(1)}\right) }} \end{aligned}$$
(82)

with

$$\begin{aligned} \mathfrak {T}_0&= {\left\langle {\ln \sum _{{x_{l}},{y_{l}}}\textrm{e}^{ \varGamma _{{{a_{l}}}{{x_{l}}}}{\left( {{\beta {{\bar{s}}}r}}\right) } +{\left( {1-\alpha }\right) }\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) } +\varGamma _{{{y_{l}}}{{b_{l}}}}{\left( {{\beta {{\bar{s}}}{{\bar{r}}}}}\right) } }}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}} \end{aligned}$$
(83)

where it is assumed that \(\beta \) is constant in \(L\). This value is then independent of \(L\) and if it is negative, the probability that \({{\varvec{b}}_L}\) is accessible from \({{\varvec{a}}_L}\) is asymptotically exponentially falling to zero. We may choose \(\alpha \in (0,1)\) as well as \(s\) and \(r\) freely except for their boundary values. But specifically for \(\alpha \) close to zero, we obtain the following expansion by matrix multiplication:

$$\begin{aligned} \mathfrak {T}_0 = {\left\langle {\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{\beta }}\right) }}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}} - \alpha {\left\langle {{\left\langle {\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) }}\right\rangle }^{s,r}_{{x_{l}},{y_{l}}}}\right\rangle }_{{a_{l}},{b_{l}}}^{p_{{}{}}} + \mathcal {O}{\left( {{\alpha ^2}}\right) }. \end{aligned}$$
(84)

At \(\beta ^*\) the zeroth order term is simply zero. The coefficient of the linear order term is exactly \(-{\mathfrak {M}^*({s},{r},{\beta })}\) and in the irregular type problem \(r\) and \(s\) can be chosen such that it is negative at \(\beta ^*\). With this choice there is then some suitable small \(\alpha > 0\), so that \(\mathfrak {T}_0\) is negative at \(\beta ^*\). \(\mathfrak {T}_0\) is continuous as a function of \(\beta \) and therefore we can then also find some \(\beta > \beta ^*\) such that \(\mathfrak {T}_0\) is still negative at the same choice of \(\alpha \), \(r\) and \(s\). This shows that the critical point \(\beta _c\) is strictly larger than \(\beta ^*\) in the irregular case, if it exists at all. \(\square \)

6 Proof of the upper side condition of Theorem 1.

In this section we derive a lower bound on accessibility, allowing us to show that the candidate threshold function in Theorem 1 does indeed satisfy the upper side of the threshold requirement for accessibility setups of regular type.

6.1 Moment bounds

To prove the lower bound on (quasi-)accessibility, we use a generalization of the second moment method. The idea of the second moment method is to bound the second moment of \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\) from above in order to apply the inequality (Alon and Spencer 2000; Hegarty and Martinsson 2014)

$$\begin{aligned} {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}> 0\right] } \ge \frac{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }^2}{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}^2\right] }}. \end{aligned}$$
(85)

\(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\) is bounded from above through the maximum length of walks and the bounded degree limiting the possible choices in each step and therefore the second moment always exists. In our proof method we do however find that, at least with our non-tight bounds on it, the second moment grows too quickly for some allele graphs to give a non-trivial bound. On the other hand, for some class of allele graphs this bound may be used to obtain a sufficient bound.

To generalize the applicability of the result, we will use a modification of the second moment method which relies on a lower order moment.

Lemma 2

Let X be a random variable over the natural numbers (including zero) with finite moment \({\mathbb {E}\left[ X^{1+\xi '}\right] }\) for some \(\xi ' > 0\), then

$$\begin{aligned} {\mathbb {P}\left[ X \ge 1\right] } \ge {\mathbb {E}\left[ X\right] }\textrm{e}^{-\frac{{\mathbb {E}\left[ X\ln {X}\right] }}{{\mathbb {E}\left[ X\right] }}} \end{aligned}$$
(86)

where an evaluation of \(0 \ln 0\) is to be taken as 0.

Proof

We know using Hölder’s inequality that for all \(0 < \xi \le \xi '\):

$$\begin{aligned} {\mathbb {E}\left[ X\right] }&= {\mathbb {E}\left[ X \mathbb {I}_{X \ge 1} \right] } \end{aligned}$$
(87)
$$\begin{aligned}&\le {\mathbb {E}\left[ X^{1+\xi }\right] }^{\frac{1}{1+\xi }} {\mathbb {E}\left[ \mathbb {I}_{X\ge 1}^{\frac{1}{1-\frac{1}{1+\xi }}}\right] }^{1-\frac{1}{1+\xi }} \end{aligned}$$
(88)
$$\begin{aligned}&= {\mathbb {E}\left[ X^{1+\xi }\right] }^{\frac{1}{1+\xi }} {\mathbb {P}\left[ X\ge 1\right] }^{\frac{\xi }{1+\xi }} \end{aligned}$$
(89)

and therefore

$$\begin{aligned} {\mathbb {P}\left[ X\ge 1\right] } \ge {\left( { \frac{{\mathbb {E}\left[ X\right] }^{1+\xi }}{{\mathbb {E}\left[ X^{1+\xi }\right] }} }\right) }^{\frac{1}{\xi }}. \end{aligned}$$
(90)

Taking the limit of \(\xi \rightarrow 0\) then completes the proof because

$$\begin{aligned} \ln {\mathbb {E}\left[ X^{1+\xi }\right] } = \ln {\mathbb {E}\left[ X\right] } + \xi \frac{{\mathbb {E}\left[ X\ln X\right] }}{\ln {\mathbb {E}\left[ X\right] }} + \mathcal {O}{\left( {{\xi ^2}}\right) } \end{aligned}$$
(91)

where the value of \(X\ln X\) is taken to be 0 by analytic continuation. \(\square \)

Because the number of walks of length \(N\) is at most exponential due to the bounded degree of \(\mathcal {A}\), while the probability of a walk to be quasi-accessible falls as fast as \(\frac{1}{N!}\), the tail of \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\) is dominated by an exponential decay. In particular all moments of \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\) exist. This allows us to apply the lemma:

$$\begin{aligned} {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\ge 1\right] } \ge {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }\textrm{e}^{-\mathfrak {K}} \end{aligned}$$
(92)

where

$$\begin{aligned} \mathfrak {K} = {\mathbb {E}\left[ \frac{\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}}{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }}\ln \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }. \end{aligned}$$
(93)

In order to prove the second part of Theorem 1 we need to show that \(\liminf {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1\right] } > 0\) with the proposed threshold function. In particular we will choose for this section \(\beta = c_L+\frac{\eta }{L}\) with some constant \(\eta \).

Per Corollary 1 then \({\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }\) converges to a non-zero value. This assures that it is sufficient to show that \(\mathfrak {K}\) does not diverge. The following method of bounding \(\mathfrak {K}\) adapts the idea used in Martinsson (2018) to account for the correlations of accessible walks using the notion of shortcuts or arcs to obtain alternative walks from a focal one.

Let \(X_\pi \) be the indicator variable that the walk \(\pi \) is quasi-accessible, then

$$\begin{aligned} \mathfrak {K} = \sum _{\pi }{\mathbb {E}\left[ \frac{X_{\pi }}{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }} \ln \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } \end{aligned}$$
(94)

where the sum is over all walks from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\) on \({\mathcal {A}^L}\) or equivalently all valid walks on \({\mathcal {A}^L}'\).

Fig. 7
figure 7

Illustration of the estimation of \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\). The short horizontal arrows (red) indicate the focal walk \(\pi \) from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\). Relative to \(\pi \) there is one non-trivial arc from \({\varvec{x}}_1\) to \({\varvec{y}}_1\), one non-trivial arc from \({\varvec{x}}_2\) to \({\varvec{y}}_2\) and two non-trivial arcs from \({\varvec{x}}_2\) to \({\varvec{y}}_3\). Between any two adjacent vertices on \(\pi \) there is furthermore one trivial arc. All other non-empty path subgraphs are not arcs, since they intersect \(\pi \) more than twice. Each walk is fully specified by choice of one of the arcs for each pair of loci. In fact it would be sufficient to choose one out-going arc per site to account for all walks. By the construction of quasi-accessibility the graph is guaranteed to be free of cycles as in the example (color figure online)

Similarly we can expand the right hand \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\) over individual walks. A graphical example for the situation described in the following can be seen in Fig. 7. We will however intentionally over-count these in the following way: Each valid walk \(\pi '\) from \(({{\varvec{a}}_L},1)\) to \(({{\varvec{b}}_L},1)\) in \({\mathcal {A}^L}'\) trivially crosses \(\pi \) in at least two vertices, namely \(({{\varvec{a}}_L},1)\) and \(({{\varvec{b}}_L},1)\). Furthermore if we list out for each valid walk \(\pi '\) the vertices it shares with \(\pi \) in \({\mathcal {A}^L}'\), then the segment of \(\pi '\) between two adjacent vertices \({\varvec{x}}\in {\mathbb {A}^L}'\) and \({\varvec{y}}\in {\mathbb {A}^L}'\) in that list does not intersect \(\pi \) a third time in \({\mathcal {A}^L}'\). We call such a segment on \({\mathcal {A}^L}'\) an arc through \({\varvec{x}}\) and \({\varvec{y}}\) on \(\pi \). An arc is said to be trivial if it is a segment of \(\pi \) itself. Immediately from the definition a trivial arc can only contain a single edge. We denote the number of non-trivial arcs which are accessible by \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }\). Each walk \(\pi '\) generates at most one arc through \({\varvec{x}}\) and \({\varvec{y}}\) on \(\pi \). Also each walk \(\pi '\) is uniquely identified by the set of arcs it generates on \(\pi \) and \(\pi '\) is accessible on \({\mathcal {A}^L}'\) if and only if all of the arcs it generates on \(\pi \) are accessible. Therefore we can bound for any valid walk \(\pi \):

$$\begin{aligned} \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \le \prod _{{\varvec{x}},{\varvec{y}}\in {\mathbb {A}^L}'}{\left( {1+\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }}\right) }. \end{aligned}$$
(95)

With this we have

$$\begin{aligned} \mathfrak {K} \le \sum _{{\varvec{x}},{\varvec{y}}\in {\mathbb {A}^L}'}\sum _{\pi }I_{{\varvec{x}}{\varvec{y}}\pi }{\mathbb {E}\left[ \frac{X_\pi }{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }} \ln {\left( {1+\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }}\right) }\right] } \end{aligned}$$
(96)

where \(I_{{\varvec{x}}{\varvec{y}}\pi }\) is an indicator variable which is 1 iff \(\pi \) contains \({\varvec{x}}\) and \({\varvec{y}}\) in \({\mathcal {A}^L}'\) and 0 otherwise.

Conditioned on the two fitness values \({F_{{\varvec{x}}}}\) and \({F_{{\varvec{y}}}}\), \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }\) becomes independent of \(X_\pi \) since \({\varvec{x}}\) and \({\varvec{y}}\) are the only vertices whose fitness values influence both the quasi-accessibility of candidate arcs and \(\pi \):

$$\begin{aligned} \mathfrak {K} \le \sum _{{\varvec{x}},{\varvec{y}}\in {\mathbb {A}^L}'}\sum _{\pi }I_{{\varvec{x}}{\varvec{y}}\pi }{\mathbb {E}\left[ {\mathbb {E}\left[ \frac{X_\pi }{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }}|{F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] } {\mathbb {E}\left[ \ln {\left( {1+\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }}\right) }|{F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] }\right] }. \end{aligned}$$
(97)

For convenience we also assume that the conditioning of the fitness values \({F_{({{\varvec{a}}_L},1)}}\) and \({F_{({{\varvec{b}}_L},1)}}\) to a difference of \(\beta \) is contained in the outer expectation.

Currently \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }\) is stochastically independent of \(X_\pi \), but still explicitly dependent on \(\pi \) in the choice of candidate arcs that need to be counted. We can remove this dependence by loosening the restriction that included arcs must not be trivial and must not intersect \(\pi \) except at \({\varvec{x}}\) and \({\varvec{y}}\). Doing so \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }\) is upper bounded by \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}}^{\prime }\), where the lack of third index indicates the loosened restriction. The resulting bound is not in general good enough for all choices of \({\varvec{x}}\) and \({\varvec{y}}\) in the sum. We will later revisit and adjust it for these special cases.

Because the logarithm is concave, the mean over it can be bounded by exchange of the two. Let \(\downarrow {\varvec{x}}\) be the projection of \({\varvec{x}}\in {\mathcal {A}^L}'\) on the first component or equivalently the 1-section of \({\mathcal {A}^L}'\). Compared to all walks on \({\mathcal {A}^L}'\) generated from walks on \({\mathcal {A}^L}\) from \(\downarrow {\varvec{x}}\) to \(\downarrow {\varvec{y}}\), arcs from \({\varvec{x}}\) to \({\varvec{y}}\) in \({\mathcal {A}^L}'\) are more restricted in the number of times vertices with projection \(\downarrow {\varvec{x}}\) or \(\downarrow {\varvec{y}}\) may or must be visited. Therefore the expectation over \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}}^{\prime }\) may be bounded by the expectation over \(\tilde{{\varvec{Z}}}_{{\downarrow {\varvec{x}}}{\downarrow {\varvec{y}}}}\).

$$\begin{aligned} \mathfrak {K} \le {\mathbb {E}\left[ \sum _{{\varvec{x}},{\varvec{y}}\in {\mathbb {A}^L}'}{\mathbb {E}\left[ \frac{\sum _{\pi }I_{{\varvec{x}}{\varvec{y}}\pi }X_\pi }{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }}|{F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] } \ln {\mathbb {E}\left[ 1+\tilde{{\varvec{Z}}}_{{\downarrow {\varvec{x}}}{\downarrow {\varvec{y}}}}|{F_{(\downarrow {\varvec{x}},1)}},{F_{(\downarrow {\varvec{y}},1)}}\right] }\right] }. \end{aligned}$$
(98)

Similarly all walks \(\pi \) through \({\varvec{x}}\) and \({\varvec{y}}\) can be separated into three segments from \(({{\varvec{a}}_L},1)\) to \({\varvec{x}}\), from \({\varvec{x}}\) to \({\varvec{y}}\) and from \({\varvec{y}}\) to \(({{\varvec{b}}_L},1)\). Each walk is uniquely determined by these three segments and for any choice of these segments forming a valid walk, their accessibility is independent under the conditioning since valid walks are selfavoiding in \({\mathcal {A}^L}'\). Taking all triples of walk segments from \({{\varvec{a}}_L}\) to \(\downarrow {\varvec{x}}\), from \(\downarrow {\varvec{x}}\) to \(\downarrow {\varvec{y}}\) and from \(\downarrow {\varvec{y}}\) to \({{\varvec{b}}_L}\), all valid walks from \({{\varvec{a}}_L}\) to \({{\varvec{b}}_L}\) through any copy of the genotypes \(\downarrow {\varvec{x}}\) and \(\downarrow {\varvec{y}}\) in \({\mathcal {A}^L}'\) are generated. This allows together with the previous arguments for the bound

$$\begin{aligned} \mathfrak {K} \le {\mathbb {E}\left[ \sum _{{\varvec{x}},{\varvec{y}}\in {\mathbb {A}^L}}\frac{ {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{\varvec{x}}}} | {F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] } {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}} | {F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] } {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{y}}}{{{\varvec{b}}_L}}} | {F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] } }{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }} \ln {\left( {1+{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}}|{F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] }}\right) }\right] }. \end{aligned}$$
(99)

Further we use that the logarithm can be bounded from above by any power law \(\ln (x) \le {\bar{\alpha }} (x-1)^\alpha \) for \(x \ge 1\), \(0< \alpha < 1\) and a constant \({\bar{\alpha }}\) depending on \(\alpha \). In particular \({\bar{\alpha }}\) as a function of \(\alpha \) can be chosen so that it is bounded except around \(\alpha = 0\), where it must diverge. Therefore, as long as we choose later any non-zero but constant \(\alpha \), the additional factor \({\bar{\alpha }}\) will not change the asymptotic order of \(\mathfrak {K}\). All in all we have the following bound, which is also represented graphically in Fig. 8:

$$\begin{aligned} \mathfrak {K} \le {\mathbb {E}\left[ {\bar{\alpha }}\sum _{{\varvec{x}},{\varvec{y}}\in {\mathbb {A}^L}}\frac{ {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{\varvec{x}}}} | {F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] } {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}} | {F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] }^{1+\alpha } {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{y}}}{{{\varvec{b}}_L}}} | {F_{{\varvec{x}}}},{F_{{\varvec{y}}}}\right] } }{{\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }}\right] }. \end{aligned}$$
(100)
Fig. 8
figure 8

Graphical representation of the bound on \(\mathfrak {K}\) given in Eq. (100). The fitness range \(\beta \) is split into three segments corresponding to paths \({{\varvec{a}}_L}\rightarrow {\varvec{x}}\), \({\varvec{x}}\rightarrow {\varvec{y}}\) and \({\varvec{y}}\rightarrow {{\varvec{b}}_L}\), indicated by straight arrows (black). Each segment contributes the given expectation value conditioned on the specified fitness difference between its endpoints into a product. The curved arrow (blue) represents the contribution of all arcs from \({\varvec{x}}\) to \({\varvec{y}}\), which contribute the given \(\alpha \)-dependent factor (color figure online)

The remaining expectation values depend only on the differences of the fitness values that they are conditioned on, not the actual placement of that difference. It is therefore convenient to use the variables \(s\) and \(r\) with \({{\bar{s}}}= 1-s\) and \({{\bar{r}}}= 1-r\) introduced previously, such that

$$\begin{aligned} {F_{{\varvec{x}}}}-{F_{{{\varvec{a}}_L}}}&= {{\bar{s}}}r\beta , \end{aligned}$$
(101)
$$\begin{aligned} {F_{{\varvec{y}}}}-{F_{{\varvec{x}}}}&= s\beta , \end{aligned}$$
(102)
$$\begin{aligned} {F_{{{\varvec{b}}_L}}}-{F_{{\varvec{y}}}}&= {{\bar{s}}}{{\bar{r}}}\beta , \end{aligned}$$
(103)

with which the outer expectation value of \(\mathfrak {K}\) can be expressed as an integral over the unit square \({\left( {s,r}\right) } \in [0,1]^2\) with a surface element \(\frac{\beta ^2{{\bar{s}}}}{2}\textrm{d}r\textrm{d}s\), which through the factor \(\frac{1}{2}\) already conditions on \({F_{{\varvec{x}}}}\) and \({F_{{\varvec{y}}}}\) being correctly ordered. We need to show that this integral is asymptotically bounded by a constant in order to show that \(\mathfrak {K}\) is asymptotically bounded by a constant from above as we intend.

The expectation value \({\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] }\) may be bounded using Eq. (72):

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}\right] } = \varTheta (1) L\textrm{e}^{L{\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }}. \end{aligned}$$
(104)

The same bound does not in general apply to the other expectation values uniformly over the integration domain due to the divergence of the constant term with vanishing fitness difference. For this reason, we split the integration region. For some sufficiently small constant \(\epsilon > 0\) we will consider integration in the regions with \(s\in [0,\epsilon ]\) and \(s\in [\epsilon , 1]\) separately and name the corresponding contributions to \(\mathfrak {K}\) accordingly with an index.

6.2 Case \(s \in [\epsilon ,1]\)

In the interval \([\epsilon ,1]\), \(s\) is bounded away from zero and therefore using Eq. (72), the expectation values \({\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}} | \beta s\right] }\) can be bounded uniformly by

$$\begin{aligned} {\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}} | \beta s\right] } \le \mathcal {O}{\left( {{1}}\right) }L\textrm{e}^{L{\varvec{\varGamma }}_{{{\varvec{x}}}{{\varvec{y}}}}{\left( {{\beta s}}\right) }}. \end{aligned}$$
(105)

For the remaining expectation values we follow the procedure used in Sect. 5 and expand with Eq. (74) to obtain a sum of locus-factorized terms. We name the special loci according to the walk segment’s index. As we have already expanded the contributions for the second segment, only the first and third remain:

$$\begin{aligned} \mathfrak {K}_{[\epsilon ,1]} = \mathcal {O}{\left( {{1}}\right) }{\bar{\alpha }}L^\alpha \sum _{l_1,l_3=1}^L{\mathbb {E}\left[ \prod _{l=1}^L\mathfrak {F}_l\right] } = \mathcal {O}{\left( {{1}}\right) }{\bar{\alpha }}L^\alpha \sum _{l_1,l_3=1}^L{\mathbb {E}\left[ \prod _{l=1}^L\mathfrak {G}_l\mathfrak {H}_l\right] }. \end{aligned}$$
(106)

Again, the usual form for loci \(l\not \in \{l_1,l_3\}\) can be given through the exponential terms in the expectation values

$$\begin{aligned} \mathfrak {G}_l = \sum _{{x_{l}},{y_{l}}\in \mathbb {A}}\textrm{e}^{\varGamma _{{{a_{l}}}{{x_{l}}}}{\left( {{\beta {{\bar{s}}}r}}\right) } + {\left( {1+\alpha }\right) }\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) } + \varGamma _{{{y_{l}}}{{b_{l}}}}{\left( {{\beta {{\bar{s}}}{{\bar{r}}}}}\right) } - \varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{\beta }}\right) }} \end{aligned}$$
(107)

and again for loci \(l\in \{l_1,l_3\}\) one or more of the exponential factors will be replaced by their derivatives. For the same reasons as used previously, in these cases \(\mathfrak {H}_l\) is uniformly bounded by a constant and therefore

$$\begin{aligned} \mathfrak {K}_{[\epsilon ,1]}&= \mathcal {O}{\left( {{1}}\right) } {\bar{\alpha }}L^{2+\alpha } \textrm{e}^{L\mathfrak {T}} \end{aligned}$$
(108)

where

$$\begin{aligned} \mathfrak {T} = {\left\langle {\ln \sum _{{x_{l}},{y_{l}}\in \mathbb {A}}\textrm{e}^{ \varGamma _{{{a_{l}}}{{x_{l}}}}{\left( {{\beta {{\bar{s}}}r}}\right) } +{\left( {1+\alpha }\right) }\varGamma _{{{x_{l}}}{{y_{l}}}}{\left( {{\beta s}}\right) } +\varGamma _{{{y_{l}}}{{b_{l}}}}{\left( {{\beta {{\bar{s}}}{{\bar{r}}}}}\right) } -\varGamma _{{{a_{l}}}{{b_{l}}}}{\left( {{\beta }}\right) } }}\right\rangle }_{l}. \end{aligned}$$
(109)

\(\mathfrak {T}\) can be considered a function of \(\beta \), \(s\), \(r\) and \(\alpha \). At \(\alpha = 0\), it is always 0 as can be verified by matrix multiplication. The first derivative towards \(\alpha \) at \(\alpha =0\) is found to be exactly \({\mathfrak {M}({s},{r},{\beta })}\). It is therefore possible to bound

$$\begin{aligned} \mathfrak {T} = \alpha {\mathfrak {M}({s},{r},{\beta })} + \mathcal {O}{\left( {{\alpha ^2}}\right) }. \end{aligned}$$
(110)

In the regular case, as \(\beta \) converges to \({\hat{\beta }}\), \({\mathfrak {M}({s},{r},{\beta })}\) is eventually bounded from above by a negative constant in the region \(s\in [\epsilon ,1-\epsilon ]\), so that in this region the integrand falls exponentially quickly to zero for suitable choice of \(\alpha > 0\), resulting in no asymptotic contribution to \(\mathfrak {K}\). In the interval \(s\in [1-\epsilon ,1]\) we need to account for the boundary term at \(s= 1\). At \(s= 1\), \(\mathfrak {T}\) is exactly \(\alpha {\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\). At the candidate threshold function \({\varvec{\varGamma }}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\) simply evaluates to \(-\frac{\ln L}{L}\) up to irrelevant higher orders in \(L\). By assumptions for the regular case we also have that the derivative \(\partial _\alpha \partial _s\mathfrak {T}\) is positive at \((s,\alpha ,\beta ) = (1,0,{\hat{\beta }})\), so that

$$\begin{aligned} \mathfrak {T} \le -\alpha \frac{\ln L}{L} + \mathcal {O}{\left( {{{\left( {\frac{\ln L}{L}}\right) }^2}}\right) } - c{{\bar{s}}}\alpha \end{aligned}$$
(111)

for some \(c > 0\). The term \(-\alpha \frac{\ln L}{L}\) exactly compensates a factor \(L^\alpha \) to the integrand of \(\mathfrak {K}\) and with a factor \({{\bar{s}}}\) in the surface element of the integration, the contribution to \(\mathfrak {K}\) from \(s\in [1-\epsilon ,1]\) is then for suitably small constant \(\alpha > 0\):

$$\begin{aligned} \mathfrak {K}_{[1-\epsilon ,1]} = \mathcal {O}{\left( {{1}}\right) }L^2\int _{1-\epsilon }^{1}\textrm{d}{s}\mathop {} {{\bar{s}}}\textrm{e}^{-cL{{\bar{s}}}\alpha } = \mathcal {O}{\left( {{1}}\right) }. \end{aligned}$$
(112)

6.3 Case \(s\in [0,\epsilon ]\)

For the integration interval \([0,\epsilon ]\) we will fix \(\alpha = 1\) and since we cannot apply the simple bound to the expectation \({\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}}|{\beta s}\right] }^{1+\alpha }\) used before uniformly in this region, we will expand it using the sum form of the expectation value. Since \(1+\alpha = 2\) now, there will effectively be two additional sums resulting from this, for which we label the corresponding locus variables \(l_{21}\) and \(l_{22}\). Again, we bring the contribution into the form

$$\begin{aligned} \mathfrak {K}_{[0,\epsilon ]} = \mathcal {O}{\left( {{1}}\right) }L^{-1}\sum _{l_1,l_{21},l_{22},l_{3}} {\mathbb {E}\left[ \mathfrak {F}_l\right] } = \mathcal {O}{\left( {{1}}\right) }L^{-1}\sum _{l_1,l_{21},l_{22},l_{3}} {\mathbb {E}\left[ \mathfrak {G}_l\mathfrak {H}_l\right] }. \end{aligned}$$
(113)

Here, since all expectation values in the numerator of \(\mathfrak {K}\) were expanded into sums, only a single factor \(L^{-1}\) remains from the expectation value in its denominator. The usual form \(\mathfrak {G}_l\) is unchanged from the region \([\epsilon ,1]\) except for the choice \(\alpha = 1\). As before \(\mathfrak {H}_l\) can be bounded by a constant for all special \(l\), but this will turn out not to be sufficient here. Suppose we used such a bound, then we would obtain

$$\begin{aligned} \mathfrak {K}_{[0,\epsilon ]} \le \mathcal {O}{\left( {{1}}\right) } L^3 {\mathbb {E}\left[ \textrm{e}^{L\mathfrak {T}}\right] } \end{aligned}$$
(114)

where \(\mathfrak {T}\) is unchanged from the previous integration region except for the choice \(\alpha = 1\). At \(s= 0\) only terms with \({x_{l}}= {y_{l}}\) can contribute to \(\mathfrak {T}\) and so it becomes 0 by matrix multiplication. The first derivative towards \(s\) can be formed directly, using that derivatives of matrix exponentials correspond to multiplication with the matrix exponent. Using that \({\left( A^{}\right) _{{{x_{l}}}{{y_{l}}}}}\mathbb {I}_{{x_{l}}{y_{l}}} = 0\) since the allele graph is simple, the derivative evaluates exactly to \(-\beta {\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\), so that:

$$\begin{aligned} \mathfrak {T} = -\beta {\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }s+ \mathcal {O}{\left( {{s^2}}\right) }. \end{aligned}$$
(115)

As \(\beta {\varvec{\varGamma }}'_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}}{\left( {{\beta }}\right) }\) is strictly positive and bounded away from zero asymptotically, this shows that \(\epsilon \) can always be chosen such that \(\mathfrak {T}\) is negative for \(s\in (0,\epsilon ]\) with negative first derivative at \(s= 0\). Consequently the integration at the boundary is of the form

$$\begin{aligned} \mathfrak {K}_{[0,\epsilon ]} \le \mathcal {O}{\left( {{1}}\right) } L^3 \int _{0}^{\epsilon }\textrm{d}{s}\mathop {} \textrm{e}^{-cLs} = \mathcal {O}{\left( {{L^2}}\right) } \end{aligned}$$
(116)

for some constant \(c > 0\). In contrast to the boundary at \(s= 1\) the surface element does not contribute here and does not yield an additional factor \(L^{-1}\). The naive bound shown above is not sufficient and two powers of \(L\) remain that we have to suppress.

To cancel these factors, we need to bound the terms \(\mathfrak {H}_l\) more carefully around small \(s\) instead of applying uniform constant bounds. In particular it would be sufficient to show that these terms introduce at least two factors \(s\) into the integrand, since the integration over \(s^n \textrm{e}^{-cLs}\) would result in a value of order \(L^{-1-n}\) instead of just \(L^{-1}\). Depending on the choices of distances between \({\varvec{x}}\) and \({\varvec{y}}\) and the choices of the special loci it is possible to provide these two factors. However not all combinations of these choices yield such a factor. The problematic cases will however turn out to be marginal in the sense that they only apply to a fraction \(\frac{1}{L}\) or \(\frac{1}{L^2}\) of the summands in the sums over special loci. Each such factor \(L^{-1}\) offsets the need for one \(s\) factor in the integrand, allowing the total contribution to \(\mathfrak {K}\) to still be constant. In the following we list all of the relevant combinations and show their contributions of \(s\) orders.

The method of bounding \(\mathfrak {H}_l\) is to consider the small-\(s\) behavior of the factors \({\left( \textrm{e}^{{s\beta }A}\right) _{{{x_{l}}}{{y_{l}}}}}\) and \({\left( A\textrm{e}^{{s\beta }A}\right) _{{{x_{l}}}{{y_{l}}}}}\) which appear in it. In particular, depending on the distance we have

$$\begin{aligned} {\left( \textrm{e}^{{s\beta }A}\right) _{{{x_{l}}}{{y_{l}}}}}&= {\left\{ \begin{array}{ll} 1 + \mathcal {O}{\left( {{s^2}}\right) } &{} {x_{l}}= {y_{l}}\\ \mathcal {O}{\left( {{s^{{d_{{{x_{l}}}{{y_{l}}}}}}}}\right) } &{} {x_{l}}\ne {y_{l}}\end{array}\right. } \end{aligned}$$
(117)
$$\begin{aligned} {\left( A\textrm{e}^{{s\beta }A}\right) _{{{x_{l}}}{{y_{l}}}}}&= {\left\{ \begin{array}{ll} \mathcal {O}{\left( {{s}}\right) } &{} {x_{l}}= {y_{l}}\\ \mathcal {O}{\left( {{s^{{d_{{{x_{l}}}{{y_{l}}}}}-1}}}\right) } &{} {x_{l}}\ne {y_{l}}\end{array}\right. } \end{aligned}$$
(118)

and due to the bounded degree of the graph, all of these bounds are uniform.

Using the bounds above, we can obtain the necessary factors of \(s\). First consider the case of all four special loci distinct. We then have for the two loci \(l_{21}\) and \(l_{22}\):

$$\begin{aligned} \mathfrak {H}_{l_{2i}}&= \frac{ \sum _{{x_{l_{2i}}},{y_{l_{2i}}}\in \mathbb {A}} {\left( \textrm{e}^{{\beta {{\bar{s}}}r}A}\right) _{{{a_{l_{2i}}}}{{x_{l_{2i}}}}}} {\left( \textrm{e}^{{\beta s}A}\right) _{{{x_{l_{2i}}}}{{y_{l_{2i}}}}}} {\left( A\textrm{e}^{{\beta s}A}\right) _{{{x_{l_{2i}}}}{{y_{l_{2i}}}}}} {\left( \textrm{e}^{{\beta {{\bar{s}}}{{\bar{r}}}}A}\right) _{{{y_{l_{2i}}}}{{b_{l_{2i}}}}}} }{ \sum _{{x_{l_{2i}}},{y_{l_{2i}}}\in \mathbb {A}} {\left( \textrm{e}^{{\beta {{\bar{s}}}r}A}\right) _{{{a_{l_{2i}}}}{{x_{l_{2i}}}}}} [{\left( \textrm{e}^{{\beta s}A}\right) _{{{x_{l_{2i}}}}{{y_{l_{2i}}}}}}]^2 {\left( \textrm{e}^{{\beta {{\bar{s}}}{{\bar{r}}}}A}\right) _{{{y_{l_{2i}}}}{{b_{l_{2i}}}}}} }. \end{aligned}$$
(119)

From Eq. (117) we can see that for all distances

$$\begin{aligned} {\left( \textrm{e}^{{\beta s}A}\right) _{{{x_{l}}}{{y_{l}}}}} {\left( A\textrm{e}^{{\beta s}A}\right) _{{{x_{l}}}{{y_{l}}}}} = \mathcal {O}{\left( {{s}}\right) } \end{aligned}$$
(120)

and therefore each of \(\mathfrak {H}_{l_{21}}\) and \(\mathfrak {H}_{l_{22}}\) contribute at least one factor \(s\), resulting in a sufficient contribution of \(s^2\) as explained above. If not all of the four loci are distinct the form of \(\mathfrak {H}_{2i}\) will be different. However, the only modifications are in the placement of derivatives of matrix exponentials. As long as still \(l_{21} \ne l_{22}\), the relevant terms which are small around \(s= 0\), namely the exponentials for the second walk segment, remain unchanged.

Therefore the remaining cases are for \(l_{21} = l_{22}\), for which we will write \(l_2\). This equality reduces the number of summands to consider by a factor \(L^{-1}\) as discussed before and consequently we need to find only one factor \(s\). In particular if either \(l_1\) or \(l_3\) are equal to \(l_2\) as well, then the weight of these cases is reduced by another factor \(L^{-1}\), so that no \(s\) is required anymore. Therefore, we can focus only on the case where \(l_1\), \(l_2\) and \(l_3\) are all distinct. For this case the contribution of locus \(l_2\) is

$$\begin{aligned} \mathfrak {H}_{l_2}&= \frac{ \sum _{{x_{l_2}},{y_{l_2}}\in \mathbb {A}} {\left( \textrm{e}^{{\beta {{\bar{s}}}r}A}\right) _{{{a_{l_2}}}{{x_{l_2}}}}} [{\left( A\textrm{e}^{{\beta s}A}\right) _{{{x_{l_2}}}{{y_{l_2}}}}}]^2 {\left( \textrm{e}^{{\beta {{\bar{s}}}{{\bar{r}}}}A}\right) _{{{y_{l_2}}}{{b_{l_2}}}}} }{ \sum _{{x_{l_2}},{y_{l_2}}\in \mathbb {A}} {\left( \textrm{e}^{{\beta {{\bar{s}}}r}A}\right) _{{{a_{l_2}}}{{x_{l_2}}}}} [{\left( \textrm{e}^{{\beta s}A}\right) _{{{x_{l_2}}}{{y_{l_2}}}}}]^2 {\left( \textrm{e}^{{\beta {{\bar{s}}}{{\bar{r}}}}A}\right) _{{{y_{l_2}}}{{b_{l_2}}}}} }. \end{aligned}$$
(121)

Following again Eq. (117), the numerator is of order \(\mathcal {O}{\left( {{s^2}}\right) }\) except if \({d_{{{x_{l_2}}}{{y_{l_2}}}}} = 1\), in which case there is a zeroth order contribution. The latter case requires additional considerations to resolve.

First, we consider the subcase with \({{\varvec{d}}_{{{\varvec{x}}}{{\varvec{y}}}}} \ge 2\). In this case it is possible that \({d_{{{x_{l_2}}}{{y_{l_2}}}}} = 1\), but if this is the case we always have another locus \(l'\) with \({d_{{{x_{l'}}}{{y_{l'}}}}} \ge 1\). Following the separation of edges in the previous section, we can handle one such locus as a special locus in exchange for another sum of order \(L\). However, the \({x_{l'}} \rightarrow {y_{l'}}\) factor contributions in \(\mathfrak {H}_{l'}\)’s numerator will then always be \(\left[ {\left( \textrm{e}^{{\beta s}A}\right) _{{{x_{l'}}}{{y_{l'}}}}}\right] ^2\) without any derivatives since \(l' \ne l_2\). From Eq. (117), such a factor results in a factor \(s^2\) compensating the additional \(L\) sum as well as the required \(s\) factor to the integrand.

Fig. 9
figure 9

Distance-1 case for \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }\). The dots represent an arbitrary acyclic subgraph. The distance-1 path indicated by the curved arrow (blue) is always accessible, assuming that the initial and final fitness values for \({\varvec{x}}\) and \({\varvec{y}}\) are correctly ordered, which we enforce through the integral boundaries. The focal path is also conditioned on being accessible. Assuming that \(\pi \) is the focal path indicated by the straight arrows (red) and \(\pi '\) the distance-1 path (blue), then in \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}\pi '}}^{\prime }\) all paths except \(\pi '\) are counted, while in \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}\pi }}^{\prime }\) at least one accessible path (\(\pi \)) is excluded, but more paths that are not arcs may also be excluded. Therefore \(\tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}\pi }}^{\prime } \le \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}\pi '}}^{\prime }\) under the stated conditioning (color figure online)

The only remaining case is then \({{\varvec{d}}_{{{\varvec{x}}}{{\varvec{y}}}}} = 1\). For this case the contribution to \(\mathfrak {K}\) is indeed not bounded as we require. However, this contribution turns out to be an overcounting issue introduced by our loosening of the arc restrictions on \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }\). Specifically, if \({{\varvec{d}}_{{{\varvec{x}}}{{\varvec{y}}}}} = 1\), we will enforce the restriction that \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime }\) should not count the direct walk segment \({\varvec{x}}\rightarrow {\varvec{y}}\) if \(\pi \) is taking this direct step. Since the direct step is always accessible given that \({F_{{\varvec{x}}}}\) and \({F_{{\varvec{y}}}}\) are ordered correctly, this segment contributes exactly 1 to the expectation value \({\mathbb {E}\left[ \tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}}}|\beta s\right] }\), which we can therefore substract from it. This is possible even if \(\pi \) does not use this direct step since \(\tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi }}^{\prime } \le \tilde{{\varvec{Z}}}_{{{\varvec{x}}}{{\varvec{y}}\pi '}}^{\prime }\) if \(\pi \) does not use the trivial arc, but \(\pi '\) does (Fig. 9). With this modification the value of \(\mathfrak {H}_{l_2}\) becomes

$$\begin{aligned} \mathfrak {H}_{l_2} = \frac{ \sum _{{x_{l_2}},{y_{l_2}}\in \mathbb {A}} {\left( \textrm{e}^{{\beta {{\bar{s}}}r}A}\right) _{{{a_{l_2}}}{{x_{l_2}}}}} {\left( A\textrm{e}^{{\beta s}A}\right) _{{{x_{l_2}}}{{y_{l_2}}}}}{\left( {{\left( A\textrm{e}^{{\beta s}A}\right) _{{{x_{l_2}}}{{y_{l_2}}}}}-1}\right) } {\left( \textrm{e}^{{\beta {{\bar{s}}}{{\bar{r}}}}A}\right) _{{{y_{l_2}}}{{b_{l_2}}}}} }{ \sum _{{x_{l_2}},{y_{l_2}}\in \mathbb {A}} {\left( \textrm{e}^{{\beta {{\bar{s}}}r}A}\right) _{{{a_{l_2}}}{{x_{l_2}}}}} [{\left( \textrm{e}^{{\beta s}A}\right) _{{{x_{l_2}}}{{y_{l_2}}}}}]^2 {\left( \textrm{e}^{{\beta {{\bar{s}}}{{\bar{r}}}}A}\right) _{{{y_{l_2}}}{{b_{l_2}}}}} }. \end{aligned}$$
(122)

Since \({\left( A^{}\right) _{{{x_{l_2}}}{{y_{l_2}}}}} = 1\), the leading order in the numerator is now \(\mathcal {O}{\left( {{s}}\right) }\), which is sufficient to obtain a bounded contribution to \(\mathfrak {K}\).

All in all, the total contributions to \(\mathfrak {K}\) are bounded in \(L\) at the candidate threshold function for the regular type, implying that there is a constant \(C > 0\), such that

$$\begin{aligned} \liminf {\mathbb {P}\left[ \tilde{{\varvec{Z}}}_{{{{\varvec{a}}_L}}{{{\varvec{b}}_L}}} \ge 1\right] } \ge C. \end{aligned}$$
(123)

This completes the proof of Theorem 1, showing that for regular setups there is indeed a threshold function \(c_L\) at which accessibility jumps in a \(\frac{1}{L}\) window from 0 to a non-zero value of at least C. \(\Box \)

It remains to improve this bound from non-zero C to \(C=1\), which we expect can be done as mentioned in Sect. 2.2, Remark 4.

7 Discussion

Cartesian power graphs provide a natural framework for describing genotype spaces composed of sequences of elements drawn from a finite set of alleles \(\mathbb {A}\). The allele graph \(\mathcal {A}\) encodes the possible mutational transitions on this set. Once the genotype-fitness map is specified according to the HoC model, which assigns fitness values to genotypes as i.i.d. continuous random variables, the rank order properties of the resulting fitness landscape are uniquely determined by \(\mathcal {A}\).

Here we have focused on the existence of fitness-monotonic paths as a measure of evolutionary accessibility (Franke et al. 2011), and proved precise results for the critical fitness difference quantile \(\beta _\text {c}\) above which accessible paths exist with positive probability. Our results quantify how accessibility increases with an increase of the number of alleles \(\mathfrak {A}\) (Wu et al. 2016; Zagorski et al. 2016), and decreases when mutational transitions are blocked or become unidirectional. For certain allele graphs, such as the path graph over three or more alleles, accessible paths do not exist for any fitness difference. Moreover, a criterion based on the behavior of Martinsson’s function \({\mathfrak {M}({s},{r},{\beta })}\) identifies allele graphs for which the behavior of the expected number of accessible paths is not informative about the existence of paths. In the words of Berestycki et al. (2017), for such allele graphs the expectation does not “tell the truth”.

The HoC accessibility problem considered here is conceptually appealing, because under the assumption of i.i.d. random fitness values, landscape accessibility is determined solely by the structure of the genotype space. However, the HoC model is not biologically realistic, as empirical fitness landscapes display varying degrees of fitness correlations (Szendro et al. 2013; de Visser and Krug 2014). The accessibility properties of correlated fitness landscapes differ significantly from those of the HoC model (Krug 2021). For example, for the much studied class of NK fitness landscapes, accessibility is determined by the structure of the interaction graph, and is low for most common structures (Hwang et al. 2018; Schmiegelt and Krug 2014). Previous work on accessibility of NK fitness landscapes has been restricted to the biallelic case, and exploring the interplay between the allele graph and the interaction graph in determining evolutionary accessibility constitutes an interesting problem for future research.