Past exam of the mathematical and Computational Biology course of the University of Cambridge / 2024

A probability density function must have integral one. For $\theta >-1$,The restriction $\theta >-1$ is essential: for $\theta \le -1$ the integral diverges and no positive normalizing constant exists.

For a moment of order $q$ with $\theta +q>-1$, integration of the probability density function givesThus the expected value and variance areThe cumulative distribution function on $(0,1)$ is $F(x)=x^{\theta +1}$. Its unique median $m$ satisfies $F(m)=1/2$, giving

For the independent and identically distributed observations, multiply the probability density functions. On the admissible parameter range $\theta >-1$, the likelihood function and log-likelihood areThese expressions presume that every observation lies strictly between zero and one; an observation outside that support would give zero likelihood under this model.

Write $s={\sum }_i\log x_i<0$. The score function isSetting it to zero gives the maximum-likelihood estimatorThis lies above $-1$. The observed Fisher information isConsequently the log-likelihood is strictly concave throughout its domain. Its stationary point is the unique global maximum, not merely a candidate local maximum. Indeed the log-likelihood tends to minus infinity as either $\theta \downarrow -1$ or $\theta \to \infty$.

The sample mean is $\bar{x}=0.7$. The method of moments equates this to the expected value from part (b):Hence the estimate is $\boxed{\tilde{\theta }=4/3\approx 1.33}$.

Here $s=\log (0.4\times 0.7\times 0.8\times 0.9)=-1.601469743$. The maximum-likelihood estimate and estimated standard error from the observed Fisher information areThe 95% Wald confidence interval is thereforeThe whole reported confidence interval lies in $\theta >-1$, though its lower endpoint is close to the boundary. With only four observations, the normal approximation underlying this Wald confidence interval need not provide accurate 95% coverage.

Keep the true parameter fixed and imagine repeatedly collecting four independent and identically distributed observations under the same model. Compute the maximum-likelihood estimate, its estimated standard error, and the Wald confidence interval afresh each time. The intended coverage statement is that approximately 95% of those random confidence intervals contain the fixed true parameter, provided the normal approximation is adequate.

After these particular observations have been collected, both endpoints are fixed. The frequentist statement does not assign 95% probability to the fixed parameter lying between these particular endpoints. Nor does the confidence interval describe the range containing 95% of future observations. The small sample makes the nominal 95% coverage an approximation whose accuracy is not guaranteed.

Both candidate models fit one parameter, but their parameter symbols have different meanings. For the exponential distribution in model B,Its second derivative is $-N/{\theta }_B^2<0$, so this is its maximum-likelihood estimate. The maximized log-likelihoods, including their normalization terms, areUse the corrected Akaike information criterion because $N=4$ is small relative to the number $k=1$ of fitted parameters:The correction is $2$ for each model, givingPrefer model A, whose corrected Akaike information criterion is smaller by $7.67252$. The ordinary Akaike information criterion also gives that ranking because the two models receive equal penalties. This compares the two proposed models; it does not establish that either is a correct description of the process. The usual corrected Akaike information criterion formula is a small-sample model-selection convention, rather than an exact finite-sample guarantee for these density families.

For the squared-denominator population map, a fixed point satisfiesThere is always the extinction fixed point $\bar{x}=0$. For a positive fixed point, divide by $\bar{x}$ and use $b+\bar{x}>0$ to obtainAt equality the candidate merges with zero and is not positive.

Differentiate the squared-denominator population map:At the positive fixed point, $b+\bar{x}=\sqrt{r}$, so its fixed-point multiplier isSince $0<b/\sqrt{r}<1$, we have $-1<F^{\prime }(\bar{x})<1$. A small perturbation ${\epsilon }_t$ obeys ${\epsilon }_{t+1}=F^{\prime }(\bar{x}){\epsilon }_t+O({\epsilon }_t^2)$, and hence shrinks locally. Thus the positive fixed point is locally asymptotically stable by fixed point stability for an iteration.

This is a stage-structured population model with a census once per maturation interval. Only adults present at census $t$ contribute to the juvenile cohort at census $t+1$; each contributes $r$ new juveniles on average. Thus $r$ is an effective fecundity, including any survival of offspring needed to reach the next census. Juveniles present at $t$ do not contribute a separate birth term during this interval.

The term $c{y}_t$ consists of juveniles that survive and mature into adults by the next census, so $c$ is their survival-to-maturity fraction. The term $(1-d)z_t$ consists of surviving adults; $d$ is the adult mortality fraction per interval. All surviving juveniles leave the juvenile compartment in one step, adults remain adults, and there is no immigration or density dependence. The reproduction term is applied to adults counted at the start of the interval, rather than only to those surviving to its end.

Eliminate the juvenile variable from the stage-structured population model. Since $y_{t+1}=2z_t$,The characteristic polynomial of this linear recurrence isso $z_t=A(3/2{)}^t+B(-1{)}^t$. The initial conditions give $z_0=0$ and $z_1=(3/4)10=15/2$, hence $A+B=0$ and $(3/2)A-B=15/2$. ThereforeThe alternating term records the initial imbalance between life stages; the growing $(3/2{)}^t$ term dominates eventually.

At a nonzero fixed point of the nonlinear stage-structured population model, the adult balance gives $d\bar{z}=c\bar{y}$. Substitute $\bar{y}=d\bar{z}/c$ into the juvenile equation and cancel $\bar{z}>0$:ThusBoth coordinates are positive precisely when $\boxed{cr>d{b}^2}$.

The Jacobian matrix at the extinction fixed point iswith characteristic polynomialUnder $cr>d{b}^2$, $p(1)=d-cr/b^2<0$, whereas $p(q)\to \infty$ as $q\to \infty$. Hence its positive eigenvalueis greater than one. Its eigenvector can be taken with both components positive: choose $(r/(b^2{q}_{+}),1{)}^T$. Thus there are biologically admissible positive perturbations that grow. The extinction fixed point is unstable by linear stability analysis for a discrete dynamical system.

Choose the diagonal matrix $W=\mathrm{diag}(w_1,\dots ,w_m)$. Multiplying a vector by $W$ scales its $i$th component by $w_i$. ThereforeThe Euclidean norm objective for ordinary linear least squares on $(WA,Wb)$ is exactly the weighted least squares objective, so the sets of minimizing vectors coincide. Zero weights simply discard the corresponding residuals.

Apply the least-squares normal equations to the transformed matrix and vector:or equivalentlyHere $W^{T}W=\mathrm{diag}(w_1^2,\dots ,w_m^2)$. A unique minimizing vector exists when the weighted design matrix has full column rank; the equations themselves remain valid without that assumption.

Using the paper's named weights $w_i$, takeThe least-squares normal equations for weighted least squares becomeMultiplying by four gives $14{\hat{\beta }}_0+33{\hat{\beta }}_1=59$ and $33{\hat{\beta }}_0+97{\hat{\beta }}_1=143$. Their determinant is $269$, givingThus, to two significant digits, $\boxed{{\hat{\beta }}_0=3.7,{\hat{\beta }}_1=0.20}$, with the slope measured in expression units per hour.

There is a wording convention to make explicit. Because the stated objective squares $w_i$, the choice above gives the final two squared residuals one-quarter of the influence of the others. If “half as much” instead means half the coefficient of each squared residual, use $W=\mathrm{diag}(1,1,1,1/\sqrt{2},1/\sqrt{2})$. Then the least-squares normal equations aregiving ${\hat{\beta }}_0=150/37\approx \boxed{4.1}$ and ${\hat{\beta }}_1=1/37\approx \boxed{0.027}$. Stating which quantity is halved resolves the ambiguity.

A BLOSUM substitution matrix is estimated from conserved, ungapped blocks of related protein sequences. Closely related sequences are clustered to reduce overrepresentation, and aligned amino acid pairs from different clusters are counted. Let $q_{ij}$ be their normalized observed pair frequency, and $e_{ij}$ the frequency expected from independent draws at the background amino acid frequencies. With unordered pairs, $e_{ii}=p_i^2$ and $e_{ij}=2p_i{p}_j$ for $i\ne j$.

For BLOSUM62, scores are rounded half-bit logarithms of these odds:A positive entry means that the aligned pair occurs more often in related sequences than the independent-background model predicts; a negative entry means less often. Adding these scores accumulates evidence across aligned positions. A score is neither a percentage identity nor a direct substitution probability.^[ref]

A diagonal substitution matrix score measures enrichment over a background match, rather than simply rewarding an identical letter. tryptophan residues are comparatively rare, so a W–W match has a small chance background frequency $p_W^2$. Its observed conservation produces a much larger ratio $q_{WW}/p_W^2$ than the corresponding ratio for the more common leucine. The BLOSUM62 diagonal therefore awards W–W more than L–L. Both the background frequency and the observed pair frequency matter; rarity alone does not determine a diagonal score.^[ref]

The number specifies the sequence-identity threshold used when clustering related sequences before deriving a BLOSUM substitution matrix. At 62%, sufficiently similar sequences are grouped and their within-cluster comparisons are excluded. Lowering the threshold merges more sequences, emphasizing comparisons across greater divergence.

Consequently BLOSUM80 is suited to closer relatives, BLOSUM45 to more divergent sequences, and BLOSUM62 is intermediate. The numbers do not specify a required identity of the query and target, nor a number of evolutionary time steps.^[ref]

Read the ten residue–residue columns from the BLOSUM62 substitution matrix. In their order in the supplied sequence alignment, their contributions areThere are three gap positions, each incurring the stated linear gap penalty of four, including the leading gap. Therefore the alignment score is $\boxed{24-3(4)=12}$.

An affine gap penalty charges a larger opening cost than extension cost, reflecting that a single insertion or deletion event can span several residues. For a gap of length $\ell \ge 1$, explicitly chooseThe supplied sequence alignment has two gap runs, of lengths one and two. Its residue–residue score remains $24$, so this choice givesMore generally, this affine gap penalty convention gives $S=24-2o-e$. There is no unique numerical rescore until opening and extension costs are chosen. A concave length-dependent gap penalty can make very long gaps progressively cheaper per added residue; distinct terminal gap penalties can accommodate incomplete sequence ends. Each choice encodes a different assumption about insertions, deletions or sequence completeness.

Use global sequence alignment, for example the Needleman-Wunsch algorithm with an appropriate substitution matrix and gap penalty. The two full-length globin proteins are related by orthology and are from closely related species, so the expected correspondence extends across both complete sequences. Unmatched ends should normally contribute to the score rather than being discarded.

With a linear gap penalty of magnitude $g$, the dynamic programming entry $F_{ij}$ is the best score for the first $i$ and $j$ residues. It satisfieswith $F_{i0}=-ig$, $F_{0j}=-jg$, and traceback from $F_{mn}$. An affine gap penalty uses separate match and gap states.

Use local sequence alignment, such as the Smith-Waterman algorithm, to compare the globin fragment with the matching region of each full-length target. A global sequence alignment would penalize long target flanks absent from the fragment. For a linear gap penalty the dynamic programming recurrence includes a restart option:Initialize boundary entries to zero, start traceback at the largest entry anywhere, and stop at zero. If the entire query is known to be valid and must be retained, fitting sequence alignment is a suitable refinement: align the whole fragment while leaving target flanks unpenalized.

First check the sequence alphabet: the prompt calls the targets “genes”. If the FASTA format file contains DNA rather than protein, translate the relevant coding sequences into amino acids before a protein-to-protein sequence alignment. If coding frames are unknown, search translated frames. A file in FASTA format specifies storage syntax, not whether the contents are nucleotides or residues.

The gain comes from avoiding a full dynamic programming table for every target. BLAST searches for short, sufficiently high-scoring word matches between query and database, then extends promising seeds. Modern versions reserve expensive gapped refinement for candidates passing earlier filters. FASTA sequence search similarly starts from short word matches, collects promising runs along alignment diagonals, and refines or joins selected regions. Word lookup and database indexing make candidate discovery inexpensive compared with considering every pair of positions.^[ref][ref]

These are heuristic algorithms: a good local sequence alignment can be missed if it has no acceptable seed, contains many interruptions, or is filtered before extension. Longer words and higher seed-score thresholds usually speed the search but reduce sensitivity; shorter words and looser thresholds recover weaker similarities at greater cost. Repetitive or compositionally biased regions can create many uninformative hits and erode both speed and specificity. No guarantee remains of finding the optimal sequence alignment for the chosen scoring scheme.

Thus the practical work depends on word length, thresholds, database composition and the number of hits. The speed-up is substantial on typical data, but does not establish a universal improved worst-case bound over $O(mn)$: an input producing many candidate hits can still require extensive refinement.

A phylogenetic epoch is a time interval between successive branching events during which the number of ancestral branches is constant. Read backward from the four sampled leaves: there are four branches until the first merger, then three until the next, then two until their common root. Write the random durations as $T_4,T_3,T_2$, with expected values $E_k=\mathbb{E}[T_k]$.

The original solution diagram below places time horizontally and labels each branch by its number of descendant leaves. Type A is the pectinate topology; type B is the balanced topology. The coloured time strips identify the three phylogenetic epochs; their drawn widths are schematic, not estimated times.

Count the number of descendant leaves on each branch during each phylogenetic epoch. From the tips backward, both topologies have branch-size lists $(1,1,1,1)$ during $T_4$, then $(2,1,1)$ during $T_3$. During $T_2$, type A has $(3,1)$ and type B has $(2,2)$.

Let $L_j$ be the total length of branches carrying $j$ descendant leaves. Adding the relevant phylogenetic epoch durations givesFor example, type A has one external branch and one three-descendant branch throughout $T_2$, while type B has two two-descendant branches. In either case the total expected branch length is $4E_4+3E_3+2E_2$, a useful accounting check. The stem above the common root is excluded because mutations there would be shared by all four samples.

In the Yule process, each of $k$ branches gives birth independently at rate $\lambda$ per generation. The next birth therefore has total rate $k\lambda$. Its waiting time has an exponential distribution, so the mean phylogenetic epoch duration isFor the ensuing calculation use $E_2=1/(2\lambda )$, $E_3=1/(3\lambda )$ and $E_4=1/(4\lambda )$. This is the inter-event waiting-time convention requested in the question, including an observation interval with four branches. Stopping the forward process immediately when the fourth branch is born would instead leave a zero terminal interval and describe a different sampling convention.

For a haploid population with $N$ possible parents, each pair of ancestral lineages merges backward at rate $1/N$ in the Kingman's coalescent approximation. With $k$ lineages there are $\left(\genfrac{}{}{0ex}{}{k}{2}\right)$ pairs. Hence the mean phylogenetic epoch duration isThus $E_2=N$, $E_3=N/3$ and $E_4=N/6$. For an autosomal sample from $N$ diploid individuals, replace the number of gene copies by $2N$, giving $4N/[k(k-1)]$. A mitochondrial sample naturally uses a haploid maternal effective population size; the common factor cancels in the ratios below.

Under the stated topology probabilities, $P(A)=2/3$ and $P(B)=1/3$. Averaging the branch lengths from part (b), the expected length producing the folded $3:1$ class isThe factor $2E_2$ in type A includes both the one-descendant and three-descendant root branches. The expected length producing the $2:2$ class isWith a common mutation rate per unit branch length, it cancels from the ratio of expected counts in this folded site frequency spectrum:For the Yule process,For Kingman's coalescent,These are ratios of ensemble expected mutation counts, not the expected value of the random ratio in one tree. The stipulated topology weights apply to both processes, and their waiting-time rates depend only on the current number of branches.

The observed folded site frequency spectrum has a larger $3:1$ to $2:2$ ratio than either expectation in part (d), and is closer to the Yule process expectation than to the constant-size Kingman's coalescent expectation. Both models can nevertheless produce deviations in one finite sample, especially because all mitochondrial sites largely share one genealogy; this observation alone is not a formal model-selection result.

For the Yule process, one possible violation is a recent reduction in the per-lineage branching rate rather than a constant $\lambda$. Lengthening the terminal four-branch interval adds one-descendant branch length and therefore preferentially increases the $3:1$ class.

For Kingman's coalescent, one possible violation is recent population expansion rather than a constant effective population size. Backward in time, coalescence is slow in the large recent population and faster in the smaller ancestral population. The resulting long external branches tend to increase the $3:1$ class relative to the $2:2$ class.

A network null model is an ensemble representing a specified baseline against which to assess an observed biological network. It should preserve the features that would otherwise confound the test: at least the number of vertices and edges, and often the degree sequence, direction, edge type, or relevant spatial constraints. It should randomize the feature whose biological organization is under investigation. An ensemble that preserves that very feature cannot test whether the feature is unusual.

For example, test whether a measured protein-protein interaction network contains more closed three-protein interaction groups than can be explained by the interaction counts of individual proteins. Define vertices as the measured proteins and an undirected edge as a reproducible detected interaction, using a fixed assay threshold and treatment of duplicates. Preselect the triangle count as the test statistic. Generate many simple random graphs by degree-preserving edge rewiring, rejecting self-loops and duplicate edges, with enough accepted swaps and separation to obtain a well-mixed reference sample. This retains the observed degree sequence while disrupting higher-order association.

Compute the same triangle count for every reference graph. Compare the observed count with this distribution; an empirical upper-tail p-value can be reported as $(1+\#\{T_{\mathrm{null}}\ge T_{\mathrm{obs}}\})/(M+1)$ for $M$ suitably sampled reference graphs. Report the size of the excess as well as its tail frequency, and check sensitivity to the assay threshold and any required connectivity constraint. A significant excess is evidence beyond the chosen degree baseline, while assay bias or spatial organization omitted from that baseline can still explain it.

A generative network model specifies a mechanism for producing networks, such as growth with preferential attachment, and asks what structures that mechanism creates. A network null model is defined by its role as a comparison baseline. These roles can overlap: the same random mechanism can be a proposed biological explanation in one analysis and a deliberately simplified null in another.

For the Erdős-Rényi model $G(n,p)$, choose the number $n$ of vertices and an edge probability $0\le p\le 1$. Visit every unordered pair of distinct vertices once and include its edge independently with probability $p$. This construction gives a simple graph with expected edge count $p\left(\genfrac{}{}{0ex}{}{n}{2}\right)$ and mean vertex degree $p(n-1)$.

Increasing $p$ at fixed $n$ makes the graph denser, increases typical vertex degree and local clustering coefficient, and tends to shorten graph distances and improve connectivity. Increasing $n$ at fixed $p$ also increases the mean vertex degree; to retain a sparse comparison one instead decreases $p$ approximately as $1/n$.

Its degree distribution is binomial distribution, approximately Poisson distribution in the sparse large-$n$ regime, so it is not a scale-free network. Connected realizations can have short average shortest path length, but their clustering coefficients are not elevated above the random baseline they themselves define. Thus they are not small-world networks under the high-clustering convention, although “small world” is sometimes used more loosely to mean only short distances.

For the Watts-Strogatz model, choose $n$ vertices on a ring, an even initial vertex degree $K<n$, and a rewiring probability $p$. Connect each vertex to its $K/2$ nearest neighbours on each side. Consider each original clockwise edge once: with probability $p$, retain its starting endpoint and replace its other endpoint uniformly, excluding self-loops and duplicate edges. The edge count remains $nK/2$, although individual vertex degrees change.

Increasing $n$ at fixed $K$ lengthens unrewired-ring graph distances. Increasing $K$ adds edges and shortens distances. Increasing $p$ introduces shortcuts, sharply reducing average shortest path length at first, then progressively reducing the clustering coefficient by destroying local triangles.

The small-world network regime occurs at intermediate $p$ with sufficient shortcuts, substantial clustering, and $K\ll n$. At $p=0$ distances are long; near $p=1$ the clustering advantage disappears. The degree distribution remains concentrated around $K$, so this is not a scale-free network.^[ref]

For the Barabási-Albert model, choose a connected seed graph on $m_0$ vertices, the number $m\le m_0$ of new edges per added vertex, and a final size $n$. A complete seed with $m_0\ge 2$ suffices. Add one vertex at a time, attaching it to $m$ distinct existing vertices selected with weights proportional to their current vertex degrees. Exclude duplicate edges and self-loops.

Increasing $n$ creates a larger network and permits older hubs to grow. Increasing $m$ makes the graph denser, with asymptotic average vertex degree $2m$, and tends to shorten graph distances. The seed affects early structure and the oldest hubs; standard linear preferential attachment gives a robust large-size tail exponent.

Its asymptotic degree distribution is $P(k)\propto k^{-3}$, so it is a scale-free network. Hubs produce short average shortest path length, but elevated clustering coefficients are not automatic: $m=1$ with a tree seed gives a triangle-free tree, and clustering decreases with size in the standard model. The strict clustered small-world network property is therefore not guaranteed.^[ref]

The supplied adjacency matrix describes a four-vertex cycle graph, with every vertex degree equal to two. Degree-preserving edge rewiring is random, so it does not determine one unique labelled output adjacency matrix. For example, replace edges $AB$ and $CD$ by $AC$ and $BD$, retaining $AD$ and $BC$. This valid switch gives, in order $A,B,C,D$,Each row still sums to two. Under the simple graph convention, all allowed outputs have the same unlabelled four-cycle topology: there are three possible labelled four-cycles. A rejected switch or later return can also leave the original adjacency matrix.

It is not a scale-free network. Its degree distribution is concentrated entirely at degree two, with no power-law tail; degree-preserving edge rewiring cannot change that.

The diagram has all possible undirected edges except $AB$: it is a complete graph on five vertices with one edge removed. Therefore $A$ and $B$ each have three neighbours, while $C,D,E$ each have four.

The neighbours of $A$ and $B$ form a triangle, so their clustering coefficients equal one. Each of $C,D,E$ has five of the six possible neighbour–neighbour edges, giving $5/6$. The nodal average shortest path length is ${\ell }_i=\frac{1}{4}{\sum }_{j\ne i}d(i,j)$. All graph distances equal one except $d(A,B)=d(B,A)=2$, henceAveraging the nodal quantities givesThe ordered-pair distance sum is $22$, and the graph diameter is two.

If the intended unnormalized “small-worldness” is the ratio of clustering to path length, the computed quantities giveThat ratio has no universal threshold for identifying a small-world network. The standard reference-normalized small-worldness, $\sigma =(\bar{C}/C_{\mathrm{ref}})/(L/L_{\mathrm{ref}})$, cannot receive a numerical value when no reference is specified and the question forbids using a network null model.^[ref]

The graph has high clustering and short paths in the descriptive sense. With nine of ten possible edges, however, those properties follow directly from its near completeness; they do not by themselves establish small-world network structure relative to a reference. This distinguishes the metrics that can be computed directly from the stronger classification that requires a stated convention.

Let ${\xi }_j$ be the displacement in step $j$, so the random walk position is $x(n)={\sum }_{j=1}^n{\xi }_j$. The single-step expected value and second moment areBy linearity of expected value, $\boxed{⟨x(n)⟩=0}$. The steps are independent random variables, so for $i\ne j$, $\mathbb{E}[{\xi }_i{\xi }_j]=\mathbb{E}[{\xi }_i]\mathbb{E}[{\xi }_j]=0$. Expanding the square gives the mean squared displacementThese are ensemble expected values; the empirical average over a finite number $N$ of independent cells fluctuates about them. The expected mean squared displacement itself does not depend on $N$.

For unbiased one-dimensional diffusion, the mean squared displacement is $⟨x^2(t)⟩=2Dt$. With $t=n\tau$, compare this to the mutant expression:Therefore the mutant diffusion coefficient is