A unifying framework for the modelling and analysis of STR DNA samples arising in forensic casework
Robert George Cowell

















	trees
	decorations.pathmorphing
	shapes,arrows

























































explExample[section]
definer[expl]Definition
opgaveExercise[section]
theoremTheorem[section]
algor[expl]Algorithm
lem[expl]Lemma
thm[expl]Theorem
prop[expl]Proposition
corr[expl]Corollary






ex
defn


alg
algno


















roman





A unifying
  framework for the modelling and analysis of STR DNA samples arising
  in forensic casework
    Robert G. Cowell
    December 30, 2023
=======================================================================================================



  This paper presents a new framework for analysing forensic DNA
  samples using probabilistic genotyping. Specifically it presents a
  mathematical framework for specifying and combining the steps in
  producing forensic casework electropherograms of short tandem repeat
  loci from DNA samples. It is applicable to both high and low
  template DNA samples, that is, samples containing either high or low
  amounts DNA. A specific model is developed within the framework, by
  way of particular modelling assumptions and approximations, and its
  interpretive power presented on examples using simulated data and
  data from a publicly available dataset.


  The framework relies heavily on the use of univariate and
  multivariate probability generating functions.  It is shown that
  these provide a succinct and elegant mathematical scaffolding to
  model the key steps in the process. A significant development in
  this paper is that of new numerical methods for accurately and
  efficiently evaluating the probability distribution of amplicons
  arising from the polymerase chain reaction process, which is
  modelled as a discrete multi-type branching process.  Source code in
  the scripting languages Python, R and Julia is provided for
  illustration of these methods.  These new developments will be of
  general interest to persons working outside the province of forensic
  DNA interpretation that this paper focuses on. 

  Keywords: DNA profiles; forensic statistics; polymerase
    chain reaction; branching process; probability generating
    functions; Fast Fourier Transform.





arabic





§ INTRODUCTION



Genetic fingerprinting, also called DNA profiling, has grown to be an
indispensable tool for identification of individuals in the
investigative and judicial process associated with criminal cases.
Judiciaries throughout the world have built up large databases of DNA
profiles of convicted, and additionally in some cases non-convicted,
individuals. These databases are based, in the main, upon
short-tandem-repeat () loci. Although newer DNA identification
techniques based upon genomic sequencing (next generation sequencing)
are being actively developed, the size of the currently existing
databases, and the comparative low cost and reliability of the
laboratory process for dealing with loci means that DNA
profiling will remain in use for several years to some.

Ever since the pioneering work instigated by <cit.>,
advances of techniques in the collection and processing of crime scene
DNA samples have led to forensic laboratories routinely processing
small amounts of DNA, for example as may be contained in just a few
cells extracted from a fingerprint. Indeed, current technology allows
for the amplification of even as little DNA as is contained within one
cell <cit.>.  Such low template DNA ()
samples can present challenges to the forensic scientist in their
interpretation. Such examples typically show signs of containing DNA
from two or more persons.  These mixtures are also subject to
various artefacts such as degradation, drop-in, drop-out and stutter
which further exacerbate the problems of their interpretation. Over
the years a variety of methods have been developed and applied to
provide interpretation in court.


It is convenient to distinguish two phases of DNA processing. The
first is the laboratory processing of the physical DNA. This part is
described in detail in <cit.>. The end
result of this is one or more electropherograms (),
described further in sec:pcr below. Each is then subject
to interpretation by the forensic statistician; interpretation issues
are described in detail in <cit.>.

In recent years probabilistic genotyping models have been
gaining acceptance within the forensic science community and in court
for interpreting challenging DNA mixtures. These models typically
describe themselves as fully continuous, by which is meant
(at least as understood by the author) that they model the peak
heights seen in the as realisations of continuous random
variables in statistical models for evaluating the weight of
  evidence of competing hypotheses presented to a court.

It is perhaps fair to say that that such models begin their
interpretation with the , or in case of the analysis of
multiple replicates from a sample, and that much of the information
contained in the processing of the DNA sample that lead to the is
not used.  It is the contention of this paper that, particularly for
low template analyses, these steps need to be explicitly or implicitly
incorporated into any statistical model for analysing DNA single
source samples and mixtures. To not do so is to lose important
information that could lead to biases or incorrect inferences. This
paper presents a framework in which such information can readily be
incorporated.

An early model for the whole process from sample to was given by
<cit.>. Their model was used only to simulate peaks, in order to relate some of the parameters of their model to the
procedures in their laboratory and observed values from runs. The
model was not concerned interpreting given . An R
package <cit.> called  <cit.>,
implements the simulation model. Recently, another model for
simulating has been proposed <cit.>.


The framework presented in this book can be thought of as an extensive
elaboration of the simulation model of <cit.>. We say
framework rather than model as many models may be
derived from the framework; in addition the framework may be extended
to cover other aspects of the processing not covered in this
paper. One key elaboration, described in more detail below, is that
the <cit.> model starts the branching process
simulation using amplicons; it ignores the fact that the processing of real DNA samples starts with genomic strands. A second
elaboration is that it does not take account of the tagging of
amplicons with dyes so that they can be observed for the . We
incorporate explicitly such the tagging in our framework. We also
extend the <cit.> model to include forward and double
backward stutters. Higher order stutters could also be included if
desired, but because of their smallness relative to other peaks in the
they are not considered important (as they would be
indistinguishable from noise).

However, not only can our framework be used to make simulated ,
but it can also be used for statistical interpretation.  A key part of
this is the use of multivariate probability generating
  functions  to model the steps in the laboratory process leading to
the . This provides a compact and elegant probability model, or
rather set of models based on approximations that are assumed. Using
standard techniques, means and variances of peak height distributions
may be found from these generating functions, and these may be used to
model peak height distributions using standard distributions, such as
normal, lognormal or gamma distributions via moment matching.  Such
distributions have been used in fully continuous models as simple and
convenient distributions: in part this has been because there is no
analytic form the for full distribution of amplicons arising in the
branching process.  However, this paper shows that such full
branching process distributions may be obtained efficiently from the
probability generating functions; there is therefore no need for
probabilistic genotyping software to assume simple standard
distributions for peak height distributions.  We show in examples that they do not
capture the intricate multi-modal distributional behaviour of
the full branching process that can occur for samples.

This paper consists of three main parts. The first introduces
background information to the objectives and problems of interpreting
forensic DNA samples.  The second part introduces elements of the
mathematical approach used in the paper, and shows how full
distributions can be efficiently computed for the branching processes
by combining a specification in terms of probability generating
functions with their evaluation using discrete Fourier transforms. The
third part presents a detailed comprehensive framework for modelling
forensic DNA problems, from which specific models may be formed by
specialisation and approximations. It describes a particular model
specialisation based on the general framework. The efficacy of the
model is exhibited with simulated and real data. 



PART: 
Background
In this part of the paper we introduce background information to the
problems we are addressing with the framework developed here.



§ DNA BACKGROUND


In this section some background information on DNA is introduced,
sufficient for the remainder of the paper; readers unfamiliar with
this background are recommended to consult
<cit.> for more
details.



 §.§ Short Tandem Repeat (STR) markers


Forensic scientists encode an individual's genetic profile using the
composition of DNA at various positions on the chromosomes.  A
specific position on a chromosome is called a locus, or
marker.  Human DNA has twenty three pairs of chromosomes:
twenty two autosomal chromosome pairs and a sex-linked pair, the X
and Y chromosomes.

The information at each (autosomal) locus consists of an unordered
pair of alleles[This is also true for the Amelogenin
  locus which occurs on both X and Y chromosomes.  However for other
  sex-linked loci the alleles might also possibly occur singly or not
  at all. For example a female does not have a Y chromosome, and so
  will not have an allele in any Y-linked locus.]  which forms the
genotype at that locus; a pair because chromosomes come in
pairs, one inherited from the father and one from the mother, and
unordered because it is not recorded from which chromosome of the pair
each allele originates.


The loci used for forensic identification have been chosen for various
reasons. Among these, we point out the two.

The first is that at each locus there is a wide variability between
individuals in the alleles that may be observed. This variability can
therefore be exploited to differentiate people.
 
 
The second reason is that, at least until recently, each (autosomal)
locus is either on a distinct chromosome, or if a pair of loci are on
the same chromosome then they are widely separated. When this occurs,
the alleles at the various loci may be treated as mutually
independent, thus simplifying the statistical analysis.  However in
recent years the numbers of loci used has increased, so that there are
now some pairs loci that are close together on the same chromosome. We
return to this issue of genetic linkage later.


The alleles of a marker are sequences of the four amino acid
nucleotides adenine, cytosine, guanine and thymine,
which we represent by the letters A, C, G and T. Each amino acid is
also called a base, and because the DNA molecule has a double
helix structure, each amino acid on one strand is linked to a
complementary amino acid on the other strand; a complementary pair of
amino acids is called a base pair.

An allele is typically named by its repeat number, usually an
integer. For example, consider the allele with repeat number 5
(commonly also referred to as allele 5 for brevity) of the marker
TH01.  This allele includes the sequence of four nucleotides AATG
repeated consecutively five times. It can be designated by the formula
[AATG]_5. Likewise allele 8 of TH01 has eight consecutive
repetitions of the AATG sequence, which may be denoted by
[AATG]_8. Repeat numbers are not always integers. For example,
allele 8.3 of TH01 has the chemical sequence [AATG]_5ATG[AATG]_3, in
which `8' refers to the eight complete four-word bases [AATG] and
the `.3' refers to the three base-long word sequence ATG in the
middle. Repeat numbers with decimal `.1' and `.2' endings are also
possible, indicating the presence of a word of one or two bases. Note
that the integer part of the repeat number counts how many complete
words of four bases make up the allele sequence but the words need not
be all identical and may vary even within loci. For example, allele 11
of the marker vWA has the base sequence TCTA[TCTG]_3[TCTA]_7. Also,
some markers are based on tri- or pentanucleotide motifs rather than
tetranucleotides as above. The base-letter sequences for many alleles
may be found in <cit.>.

When the repeat numbers of the two alleles of an individual at a
marker are the same, then the genotype for that marker is said to be
homozygous; when the repeat numbers differ, the genotype for
that marker is said to be heterozygous.

The repetitive structure in the alleles gives rise to the term
short tandem repeat () marker to describe these loci; they
also go by the name of microsatellites.  Note that for other
purposes of genetic analysis it is common to use single-nucleotide
polymorphisms (SNPs) which are defined as DNA sequence variations that
occur when a single nucleotide (A, T, C, or G) in the genome sequence
is altered.


Forensic identification using s is based upon the size of the
allele. However, for some loci, different people might not have the
same nucleotide sequence for one or more alleles of a specific repeat
number. This information is not used by the framework developed in
this paper, but is the basis of finer genetic discrimination available
to the next-generation sequencing methods currently being developed
that were mentioned in the introduction.


Within a population the various alleles of markers do not occur
equally often, some can be quite common and some quite rare.  When
carrying out probability calculations based on DNA, forensic
scientists use estimates of probabilities based on allele frequencies
in profiles of a sample of individuals. The sample sizes typically
range from a few hundred to thousands of individuals.  For example,
<cit.> presents tables of US-population -allele
frequencies for Caucasians, African-Americans, and Hispanics based on
sample sizes of 302, 258 and 140 individuals. These are used to
estimate the genetic profile probabilities of an individual: for the
autosomal loci and the locus Amelogenin, independence of loci
expressed through the lack of genetic linkage means that they may be
found by multiplying probabilities across the loci.  Although the process described below applies to other sex-linked loci, the linkage
of loci on the sex-linked loci presents special problems in estimating
genetic profile probabilities. This issue will be addressed in another
paper, although some of the discussion below will not be specific to
autosomal loci.




 §.§ The PCR process



The DNA collected from a crime scene for forensic analysis consists of
a number of human cells from one or more individuals. Note that each
cell of an individual will contain two alleles (diploid cells) for
each autosomal marker, whereas sperm cells have only one allele
(haploid cells).  This means that in a mixture, a particular
individual will contribute the same number of alleles for each
marker. In order to identify the alleles that are present, a DNA
sample is first subjected to chemical reagents that break down the
cell walls so that the individual chromosomes are released into a
solution. A small amount of this solution is used to quantify the
concentration of DNA; the typical unit of measurement is picograms per
microlitre, the DNA in a single human cell having a mass of between
six and seven picograms. Having determined the density of DNA in the
sample, a volume is extracted that is estimated to contain a certain
quantity of DNA, typically around 0.5 nanograms, equivalent to around
75 human cells.  The DNA in the extract is then amplified using the
polymerase chain reaction (PCR) process. This involves adding
primers and other biochemicals to the extract, and then
subjecting it to a number of rapid heating and cooling cycles. Heating
the extract has the effect of splitting apart the two complementary
strands of DNA, the cooling phase then allows free floating primers
and amino acids to bind with these individual strands in such a way
that the DNA is copied. By the action of repeated heating and cooling
cycles, typically around 28 altogether, an initially small amount of
DNA is amplified to an amount large enough for
quantification. Mathematically, the amplification process may be
modelled as a branching process
<cit.>. The amplification
process is not 100% efficient, that is, not every allele gets copied
in each cycle. This means that if two distinct alleles in a marker are
present in the extract in the same amount prior to amplification, they
will typically occur in different amounts at the end of the PCR
process.


Note that after breaking down the cell walls to release the genomes,
there could be sufficient DNA in the sample to carry out several
independent PCR amplifications with sub-extracts. When this is done it
is called a replicate run.


To understand the quantification stage of the post PCR amplified DNA,
it is important to know that the amplification process does not copy
only the repeated DNA word segment of a marker, it also copies DNA at
either end. These are called flanking regions, and their
presence is important in performing the PCR process. Regions at either
end of the flanking region, called primer binding regions,
are where the primers bind to the DNA to initiate copying. These
regions are labelled as 5' and 3'; there is one for each of the two
strands making up the double helix of the DNA. During the heating
cycle the double helix separates into two strands, which for reasons
that will become apparent later we denote by the two genomic strands
by g and g_d. During the cooling phase one of the strands primers
will bind at the 5' end, and on the other at the 3' end. After the
primers have attached the Taq polymerase adds individual bases to
complete the copying. This leads to two strands that we call h and
h_d, (with h for half-strand).

The process is illustrated in the following sequence of figures. In
the first figure we have a fragment of DNA consisting of two
complimentary joined strands g and g_d which are long enough to
contain the repeat structure of an allele of interest and primer
binding regions on either side of the allele. The `teeth' in the
figure indicates the individual bases making up each strand.



  
    < g r a p h i c s >




Heating breaks the bonds of the base-pairs of the DNA molecule, so
separating the strands.  During cooling primers,indicated by the
smaller sets of teeth. bind to each of the 5' ends of the strands.
Note that the sequences of bases in the primers attaching to the two
ends are in general different, because of the differences in the base
sequence in the primer-binding regions of the two strands.



  
    < g r a p h i c s >




After the primer binds at one end the Taq polymerase extends to copy
to beyond the primer region at the other end.  The g strand
generates the h_d strand for its copy, and the g_d strand the h
strand for its copy.



  
    < g r a p h i c s >




During the next heating cycle the g/h_d strand pair separate, and
the g_d/h pair separate. The g and g_d strands repeat the
process as above. During the cooling phase, primers bind to the 5'
ends of the h and h_d strands,



  
    < g r a p h i c s >




and the Taq polymerase adds bases to the ends, thus making with the
h strand a complementary a_d strand, and with the h_d strand, a
complementary a strand.


  
    < g r a p h i c s >





During the next heating cycle the h/a_d strand pair separate, as
does the h_d/a strand pair.  The h and h_d strands behave as
above (as the g and g_d strands continue to do so as
well). Primers bind at the 5' ends of the a and a_d strands,




  
    < g r a p h i c s >




and the Taq polymerase completes the copy, so that the a strand
makes a complimentary a_d strand, and the a_d strand makes a
complimentary a strand.


  
    < g r a p h i c s >





Thus an amplified allele will consist of the allele word repeat
sequence region and two flanking regions, and will have a length
associated with it which is measured in the total number of base pairs
included in the word sequence and the flanking regions.  This is
called an amplicon.


  
    < g r a p h i c s >




Further heating cycles lead to an exponential growth in the number of
a and a_d strands with the number of cycles. After around 28 or 29
cycles there is a sufficient number of amplicons to measure their
amount.


For each marker, the DNA sequence (hence the size) forming each of the
two flanking regions is constant, but different across markers. Thus
quantifying a certain allele is equivalent to measuring how much DNA
is present of a certain size. This is carried out by the process of
electrophoresis, as follows.

During the binding of primers a fluorescent dye is also attached. We
said earlier that the primers for the two ends of the strand are
different. Typically dyes are attached to only of  the two types or
primers: we can now reveal that the d subscript used in the
description indicates the presence of a dye tag, thus a_d represents
the strand in the complementary pair of strands making up the
dye-tagged amplicon, and a other strand that is not tagged with a
dye. (Note that h_d represents a dye-tagged half-strand, but that
g_d strands are not dye tagged. However it simplifies notation to
treat the g_d as if they are, rather than introduce another symbol
to distinguish the g strand from its compliment.)  For brevity in
the remainder of the paper, we shall refer to a_d as a dye-tagged
amplicon, and a and an untagged amplicon.

Hence we wish to measure the number of dye-tagged amplicons,
a_d. The description above has assumed 100% efficiency in the
amplification process. However, not every strand will get a primer
attached in each cycle, and the process will operate at less than
100% efficiency. Typical amplification efficiencies are in the range
0.8-0.9.


Several colours of fluorescent dyes are used to distinguish similarly
sized alleles from different markers.  The amplified DNA is drawn up
electro-statically through a fine capillary to pass through a light
detector, which illuminates the DNA with a laser and measures the
amount of fluorescence generated. The latter is then an indication of
the number of alleles tagged with the fluorescent marker. The longer
alleles are drawn up more slowly than the shorter alleles, however
alleles of the same length are drawn up together. This means that the
intensity of the detected fluorescence will sharply peak as a group of
alleles of the same length passes the light detector, and the value of
the intensity will be a measure of the number of alleles that
pass. The detecting apparatus thus measures a time series of
fluorescent intensity, but it converts the time variable into an
equivalent base pair length variable. The data may be presented to a
forensic scientist as an electropherogram () as illustrated
by the simulated shown in fig:epg. Each panel in the EPG
corresponding to a different dye. The horizontal axes indicate the
base pair length, and the vertical axis the intensity.


In the absence of artefacts, a peak in the indicates presence of
an allele in the sample before amplification. The peak height is a
measure of the amount of the allele in the amplified sample expressed
in relative fluorescence units (RFU). The area of the peak is
another measure of the amount, but this is highly correlated with the
height <cit.>. Both peak height and peak area are
determined by software in the detecting apparatus.

We shall call the peak size information extracted from the the
profile of the DNA sample, or more briefly, the DNA
  profile. Commonly, DNA profile also refers to the combined genotype
of a person across all markers.


In measuring the peak heights, low level noise give rise to small
spurious peaks. A peak amplitude threshold, called the
analytic threshold, may be set by the forensic analyst
whereby peaks below the analytic threshold level are ignored.  Thus an
allele present in the DNA sample will not be recorded as observed if
the peak it generates is below the analytic threshold; when this
happens a dropout of the allele is said to have occurred. A
dropout of an allele can also occur if no genomic strands containing
the allele are selected for amplification.  Dropout is an artefact
that can make the analysis of DNA samples difficult. Another common
artefact is stutter, where­by an allele that is present in
the sample is mis-copied at some stage in the PCR amplification
process, and (for a tetraneucleotide marker) a four base pair word
segment is omitted.  This damaged copy itself takes part in the
amplification process, and so yields a peak located four base
pairs[For tetrameric loci.]  below the allele from which it
arose.  More rarely, two repeats are omitted during a cycle,
which is called double stutter, or an extra repeat is
inserted, which is called forward stutter; once formed these
artefacts can themselves replicate in subsequent cycles.

Another artefact is known as dropin, referring to the
occurrence of small unexpected peaks in the EPG. This can for example
be due to sporadic contamination of a sample either at source or in
the forensic laboratory.


Finally, a mutation in the flanking region can result in the allele
not being picked up at all by the PCR process, in which case we say
that the allele is silent. An allele can also be undetectable
and thus de facto silent because its length is off-scale and
the peak does therefore not appear in the . Note that an allele
might be silent for a kit made by one manufacturer but not another;
this is because different manufacturers of kits for performing use different primer binding regions.




§ FROM SAMPLE TO EPG


In this section we give a description of the process of going from a
DNA sample, recovered for example from a crime scene, to the
electropherogram () that is used to make inferences about the
constitution of the DNA in the sample. For a much more detailed
account of the laboratory processes please see
<cit.>.

We start with a descriptive summary, and follow this with a more
detailed mathematical specification.



 §.§ Steps in the process

The following is a summary of the sequence of steps taken to obtain an
from a DNA sample, recovered for example at a crime scene;
fig:epgsteps illustrates the sequence of steps involved.




  * Collect sample of DNA. The sample might be a single source
  trace, that is, it contains DNA from a single individual, or it
  could be a mixture, that is a trace with DNA originating from two or
  more people.

  * Take a sub-sample of the sample and using bi-chemical reagents
  extract nuclear DNA to create an aliquot (solution of extracted DNA-
  check this is the correct meaning). Typically this will be in a
  mini-tube of a volume of between 50 to 100, depending on
  laboratory standard operating procedures.

  * Take a small extract of the aliquot, typically around 2,
  to determine the concentration of DNA in the aliquot. Usually this
  will be carried out using qPCR
  <cit.>, so that (if possible) an optimal
  amount DNA can be used for the step.


  * Based upon the estimated density of DNA in the aliquot, take
  enough of the aliquot for amplification so that there is an
  optimal amount of genomic DNA in the mini-tube for the amplification is to be carried. For samples, obtaining a
  sufficient amount might not be possible.  Typically a maximum volume
  of 20 is taken from the aliquot, so that there is sufficient
  volume left in the PCR mini-tubes for the PCR primers. In some cases
  the aliquot is split into several parts called replicates
  and each is separately amplified via .
  

  * Carry out the amplification of each replicate using some
  kit and protocol.  Usually, but not always, the same manufacturer
  kit is used for all of the replicates when more than one replicate
  is made.


  * An is then produced from each replicate by capillary
  electrophoresis.










§ THE OBJECTIVES OF AN ANALYSIS


In the previous sections we have given a brief qualitative description
of how we obtain an from a DNA sample.  In this section we
describe the information a forensic scientist wishes to obtain from an
. In later sections we develop a mathematical framework to help
realise those aims.

A forensic scientist may wish to use the information in an to
calculate the weight of evidence for two competing hypotheses
regarding the contributors to the DNA sample; another objective is to
try and identify the likely genotypes of contributors
<cit.>, a
process known as deconvolution.

We consider each of these objectives in turn.


  
Weight of evidence
 

In an adversarial court setting, we have two competing hypotheses, one
argued for by the prosecution, called the prosecution
  hypothesis _p, the other called the defence hypothesis
_d. Note that the defence is not obliged to propose an
alternative to the prosecution hypothesis _p; usually the
defence hypothesis is proposed by the prosecution for comparative
purposes.  The available evidence E in the case consists of
the peak heights as observed in the as well as the set of
genotypes of some known individuals; denote the known genotype
information by G. The prosecution and defence hypotheses differ in
their assumptions as to whose DNA is in the sample that produced the
. In most cases where a suspect S is on trial, the prosecution
case would be that S contributed to the sample and the defence case
that S did not. However this is not always the case. For example, a
DNA sample might be recovered from a swab of an area of skin of the
suspect where physical contact is alleged to have been made during an
assault on a victim V. Under this scenario, the presence of S's
DNA in the sample is not disputed, it is the presence of the victim's
DNA offered by the prosecution as proof of the assault that is
disputed.


The strength of the evidence <cit.>
is normally represented by the likelihood ratio:

    LR=L(_p)/L(_d)=( E _p)/(E _d).



This may be expressed on a base-10 log scale of units called the
ban, introduced by Alan Turing <cit.> so that one ban
represents a factor 10 on the likelihood ratio.  Then the weight
  of evidence is =log_10LR in bans
<cit.>.



The numerator and denominator in the likelihood ratio are calculated
based on models having the form

    ( E , G)= ∑_𝐠( E 𝐠) ( 𝐠, G),

so that the model for the conditional distribution
( E 𝐠) of the evidence given the genotypes
𝐠 of all contributors is the same for both hypotheses,
whereas the hypotheses differ concerning the distribution
( 𝐠, G) of genotypes of the contributors. Note
that the genotypes 𝐠 that are summed over may include those
of one or more unknown contributors, that is, individuals
whose genetic profiles are not known (not included in G).

Note that ( E 𝐠) has a dependence, not shown, on the
amount of DNA from each contributor in the sample, and other factors
such as sample degradation: this implies a similar implicit dependence
for the left hand side.

 


  
Deconvolution of DNA mixtures

In the deconvolution of a DNA sample, typically a mixture, we assume
that there are one or more genetically untyped contributors to the DNA
sample. We then wish to find the genotypes of these untyped
individuals. Typically how this is done is that for each individual
separately a ranked list of genetic profiles ordered by their
likelihood is produced.  The potential profiles having high likelihood
could, for example, be checked against an offender database for a
match.  However, sometimes a ranking of the joint genotypes of two or
more individuals may also be of interest.



PART: 
Mathematical formulation

In this part of the paper we give the mathematical notation
underpinning the framework of this paper.  From the right hand side of
(<ref>) we see that there are two components of the
framework that require specification.  One part is specifying the
(conditional) probability for the data, ( E 𝐠),
the other specifying the genetic profile probabilities,
( 𝐠, G).  As the latter is relatively
uncontroversial, and is common to all models for the
( E 𝐠) distribution, we discuss that first. We then
discuss the simulation model of <cit.>, and show how
the full distribution may be obtained without simulation.  We then
extend the analysis to finding the distribution for the more realistic
model in which we start with genomic strands and we find the
distribution of dye-tagged amplicons. We then elaborate this in stages
leading to the general framework that includes background noise,
drop-in, forward and double-reverse stutters, degradation and
inhibition, for possibly multiple replicates analysed with possibly
multiple kits from possibly multiple independent samples.




§ SPECIFYING GENETIC PROFILE PROBABILITIES


Given the genotypes of contributors, the framework formulated in this
paper for specifying peak-height likelihoods is applicable to any type
of STR locus, that is for evaluating ( E 𝐠). However, if there are untyped contributors, evaluation
of ( 𝐠, G) can be problematic for sex-linked
loci. In this paper we shall assume that the set of loci are all
autosomal loci plus, optionally, Amelogenin.  We shall also assume
that all contributors are unrelated. With these assumptions
(<ref>) simplifies to a product over the loci L, with
𝐠_𝐥 denoting a genotype in the set G_l of genotypes on
the locus l∈ L:


    ( E , G)= ∏_l ∈ L(∑_𝐠_𝐥( E 𝐠) ( 𝐠_𝐥, G_l)),


We assume that the population from which each contributor comes from
is known, and that allele frequency estimates are available. We do not
assume that the contributors all come from the same population. It is
assumed that populations substructure correction parameters are known
for each population. Optionally, the finite size database correction
of <cit.> may be applied. We do not have to be
concerned with linkage between autosomal loci because of the
assumption that contributors are unrelated. Relatedness amongst
contributors could be incorporated, but at the cost of complicating
the presentation presented here, especially if there are linked loci.
For a detailed discussion of these issues relating to populations and
genotype probability estimation, as they relate to forensic
applications, see
<cit.>.
A companion paper by the author is planned that proposes solutions to
the twin problems of Y-STR haplotype probability estimation and the
resolution of Y-STR mixtures.





§ THE SIMULATION MODEL OF <CIT.>: THE AMPLICON
  MODEL


In this section we summarise the simulation model of
<cit.>, for further details please see the original
paper. We shall refer to this model as the amplicon model.

The model presented in <cit.> has variants for diploid
cells and haploid cells. As mentioned earlier, their model fails to
take account of the fact that the process starts with genomic material
rather than amplicons, and it also ignores the dye-tagging of the
amplicons during the process. Despite these short-comings, the
model provides a convenient starting point for the more realistic
models presented later.  We start off with the case of no stutters;


 §.§ Simulating the process without stutters


The model of <cit.> is a simulation model. It assumes
that initially we have a number n_c of cells in our small subsample
of fig:epgsteps.  We concentrate on just one allele from one
locus. If the cell is diploid, and the individual is homozygous, then
there will be 2n_c such amplicons within the cells to start with;
otherwise if the individual is heterozygous on the locus, or the cell
is haploid, then there will be n_c such amplicons in the cells to
start with.  For whichever case holds, denote by N the total number
of amplicons initially.

In the first stage the amplicons are extracted from the cells. This
process is not 100% efficient, more typically only between 10-30%
are extracted intact <cit.>. Following
<cit.>, let π_eff denote the extraction
probability that an individual amplicon is released intact into the
aliquot. Assuming independence of the release of the distinct
amplicons, the total number of intact amplicons in the aliquot is
therefore binomially distributed with distribution
(N,π_eff).

A fraction of the aliquot π_aliquot is then taken for .  If
we let ϕ = π_effπ_aliquot, then the total number of
amplicons intact and selected for amplification will be Binomially
distributed as
(N,π_effπ_aliquot)≡(N,ϕ).
 
These are now subject to amplification, where it is assume that
in each cycle, each amplicon makes a copy of itself with probability
π_pcr If we assume that there are K cycles, then the total
number of amplicons can be simulated using the following algorithm,
starting with N amplicons:
 

[Simple amplicon simulation]
  
  
  
  * Randomly sample n from (N,ϕ).
  
  * For K times do:
    
    
  * Sample m from (n,π_pcr)
    
  * Update n := n+m
    
  
  

This algorithm is very simple to implement. If it is run many times, a
histogram, or kernel density estimate, plot may be made of the
distribution of the number of amplicons.


Python code presented in sec:py:ampkd produced the kernel
density plot shown in fig:ampkd, based on a million
simulations.  The parameter values used were π_eff = 0.6,
π_aliquot = 20/66, π_PCR = 0.8 and K=28 cycles (values
taken from the <cit.> paper), for a single starting
amplicon, N=1.  Note that for each simulation there is a probability
ϕ = 0.6× 20/66 = 2/11 that the amplicon is selected for
amplification , so that we have a total dropout probability of
9/11. As this will dwarf the rest of the plot, these zeros have been
removed before the kernel density has been estimated.







 §.§ Generating the full distribution (without stutters)


We now show how the simulation model can be represented by using
probability generating functions, following the approach of
<cit.> based on Galton-Watson cascade processes.  Consider
an initial single amplicon.  In the first PCR cycle it can either
amplify with probability π_PCR, or fail to amplify with
probability 1-π_PCR. The probability distribution of the number
of alleles after one cycle can therefore be represented by the
probability generating function ():


    F_1(t) = f(t) = (1-π_PCR)t + π_PCR t^2.

After two cycles the number of alleles has the 
    F_2(t) = f(f(t)),

and after 3 cycles it has the 
    F_3(t) = f(f(f(t))).


More generally, we have


  
  The for the number of molecules after r cycles, given
  that there is exactly one prior to any amplification cycle, is
  
    F_r(t) = f(f(f(… f(t) … )))

  where
  
    F_1(t) = f(t),    F_s+1(t) = f(F_s(t)),  (s=1,2,3,…).



This is Theorem 1 of Good <cit.> (but here specialised to
F(t) given above).
 
The above gives the for a single starting amplicon: if there are
K cycles the will be a polynomial in t of degree
2^K. However, for the simulation model we start not with a amplicon,
but with N initial amplicons, that are each sampled independently
for amplification with probability ϕ.  The taking into
account the initial number and the pre-sampling is
<cit.>

    (1 - ϕ + ϕ F_K(t))^M

a polynomial of degree M× 2^K. This is simply the functional
composition of the Binomial (1-ϕ + ϕ t)^N for the
pre-sampling of the amplicons with the of the number of
amplicons arising from a single amplicon in the branching
process.
 
Let us rewrite the recursion in (<ref>) on substituting
f(t); for later convenience we also replace π_PCR by p_t,
thus the recurrence relation becomes:


    F_s+1(t) = (1- p_t)F_s(t) + p_t F_s(t)^2.

 
with initial value F_0(t) = t.  It is in principle possible to find
the polynomial F_K(t) quite simply by iteration using a computer
algebra package, however the growth in the numbers of terms in the
polynomials means that there is a relatively small limit to the number
of iterations that can be carried out before computer memory is
exhausted of the order of K=10 or so, even if numerical values for
p_t are used. In addition in the later stages, the brute force
evaluation of the quadratic term F_s(t)^2 grows quadratically in
complexity, thus making the later stages of iterations slower and
slower.


There is however another way to evaluate the numerically, by
noting that the quadratic expression F_s(t)^2 is simply the of
the convolution of two (identical) probability distributions. The
convolution can be carried out using a Discrete Fourier
  Transform (), which may be done efficiently using a
Fast Fourier Transform (). Moreover, as we shall show,
the binomial sampling composition can also be carried out using the
.  Before given the details, the reader may care to look at
fig:ampfft28, which shows the exact distribution for a single
starting amplicon evaluated using the .







  §.§.§ The discrete Fourier transform

Before presenting the algorithm for generating the full distribution
using the , we first give a brief review of the .  Suppose
that a sequence of (real) numbers (x_0, x_1, x_2 , …, x_n-1) is
given, and an integer N ≥ n.  Then the sequence of numbers
defined by

    X_k = ∑_j=0^N-1 x_je^-2π i jk/N,

where i = √(-1), is the discrete Fourier transform of the
sequence (x_0, x_1, x_2 , …, x_n-1). There is an inversion
formula:

    x_k = 1/N∑_j=0^N-1 X_je^2π ijk/N


The key result that we use is that multiplication of two large
polynomials may be done efficiently using the , which may be
calculated efficiently using the . Suppose that x(t) is a
polynomial of degree n in t with coefficients
(x_0, x_1, x_2 , …, x_n), and that y(t) is a polynomial of
degree m in t with coefficients (y_0, y_1, y_2 , …,
y_m). Let z(t) = x(t)y(t); this is a polynomial of degree m+n
with coefficients (z_0, z_1, … z_n+m) where

    z_k = ∑_j=0^k x_jy_k-j,

(defining x_j=0 for j > n and y_k-j=0 for k-j >m).  Let N
be an integer such that N ≥ n+m+1. if we extend the sequence
(x_0, x_1, x_2 , …, x_n) with zeros to create a new sequence
(x_0, x_1, x_2 , …, x_N-1), and similarly extend the sequence
(y_0, y_1, y_2 , …, y_m) with zeros to create a new sequence
(y_0, y_1, y_2 , …, y_N), then we may form the of these
extended sequences:

    X_k    = ∑_j=0^N-1 x_je^-2π i jk/N
    
      Y_k    = ∑_j=0^N-1 y_je^-2π i jk/N

If the sequence (z_0, z_1, … z_n+m), extended with zeros if
required to make the sequence (z_0, z_1, … z_N-1), then the
of this sequence is

    Z_k = ∑_j=0^N-1 z_je^-2π i jk/N

and we have that

    Z_k = X_k Y_k    k ∈{0,N-1}.

so that

    z_k = 1/N∑_j=0^N-1 X_jY_j e^2π ijk/N

Hence to multiple the polynomials x(t) and y(t), we choose N
sufficiently large, take the s of the coefficients of x(t) and
of y(t), multiply the two s term-wise, and then take the
inverse . The computation of the forward and backward may be
carried out efficiently using the algorithm, an algorithm by
Gauss that dates back to 1805 (see <cit.> for an
interesting history of the , and <cit.> for a recent
monograph on algorithms.).


We can now present an algorithm for generating the full distribution
using the . It is given in alg-gillampliconfft for the
case of a single starting amplicon and no binomial presampling.  Let
F[] be a vector with indices starting from zero, such that F[n]
denotes the coefficient of t^n in the . With K cycles, the
number of amplicons will range up to 2^K, hence F[] must be a
vector of size at least 2^K+1.


[Single amplicon distribution using the ]
  
  
  
  * Initialise F[] to be a vector of size at least 2^K+1,
    with all entries 0 except F[1]=1.
  
  * Let F[] denote the of F[] .

  
  * For K times do:
    
    
  * for each element F[f] of F[] update
      F[f] := (1-p_t) F[f]+ p_t  F[f]^2
    
  
  * Set F[] equal to the inverse of F[].
  
  

Note that taking the is a linear operator; denote it by L and
its inverse by L^-1.  Hence we may write (with × denoting
element wise multiplication)

    F_s+1(t)    =    (1- p_t)F_s(t) + p_t F_s(t)^2
       =    (1- p_t)F_s(t) + p_t L^-1( L(F_s(t)) ×  L(F_s(t)) )
       =    (1- p_t)L^-1(L(F_s(t)) + p_t L^-1( L(F_s(t)) ×  L(F_s(t)) )
       =    L^-1((1- p_t)L(F_s(t)) + p_t  L(F_s(t)) ×  L(F_s(t)) ),

thus justifying alg-gillampliconfft.

Taking account of starting with M amplicons, and binomially sampling
them with probability ϕ is almost trivial. We have, using the
linearity of L and its inverse:

    (1 - ϕ + ϕ F_K(t))^M    =    L^-1 L( (1 - ϕ + ϕ F_K(t))^M) )
       =   L^-1( (L(1 - ϕ )+ ϕ L(F_K(t))^× M)

where the superscript ^× M denotes taking the M-th power of
each element in the transform. Hence to take this into account we
modify alg-gillampliconfft thus, (extending the initial size
of F[] for the higher number of amplicons that could result):

[Amplicon distribution with binomial pre-sampling
    using the ]
  
  
  
  * Initialize F[] to be a vector of size at least M(2^K+1),
    with all entries 0 except F[1]=1.
  
  * Let F[] denote the of F[] .

  
  * For K times do:
    
    
  * for each element F[f] of F[] update
      F[f]:= (1-p_t) F[f]+ p_t  F[f]^2
    
  
  * for each element F[f] of F[] update
    F[f]:= (1-ϕ + ϕ F[f])^M.
  
  * Set F[] equal to the inverse of F[].
  
  

Now the extreme right-hand tail of the amplicon distribution has very
low probabilities, and for all practical intents and purposes can be
taken to be zero. Hence we may take the F[] vector to have size
M 2^K which makes finding the using the much more
computationally efficient, especially if M itself is a power of
2. Using this approximation, alg-gillampliconfftbin is
readily implemented.  r-alg-gillampliconfRcode presents the
seven(!) lines of code R code for this, repeated here.






fig:targetmargK16M50Phi06 shows the amplicon distribution
obtained using the code given above for a low template amount, in
which multimodality is clearly present, and would be missed by current
probabilistic genotyping software based assuming a unimodal
probability distribution (such as lognormal or gamma) for peaks
heights.[Equivalent Python code, used for generating
  fig:ampfft28 is given in sec:py:ampfft. Running
  the Python code for the full 28 cycles took around 3 minutes on a
  laptop with an Intel i7TM processor, and used
  around 30Gb of ram.]


It is also worth emphasising that the model naturally includes
dropout, and that the full distribution calculated using
alg-gillampliconfftbin enables the means to assess dropout
probabilities –there is no need to posit a separate additional model
for allelic dropout, such as the logistic regression model (or
variants thereof) of <cit.>.  Indeed, in a
follow-up paper, <cit.> refine their logistic
regression approach with a probit model based on and compared with the
amplicon model of <cit.>.


To illustrate this, we consider an amplification probability of
p=0.85. To keep computations manageable we take the number of cycles
to be K=15, take ϕ = 0.1, and set the analytic threshold to be
40000 amplicons. fig:dropoutamp shows the dropout
probabilities P(D) as a function of the number of starting cells -
note that the horizontal axis is on a log-scale for comparison with
Figure 1 of <cit.>. The plot shows two curves,
with the red corresponding to a homozygous individual, and the black
to a heterozygous individual. (R code to generate this plot is given
in r-fig:dropoutamp.)







 §.§ Alternative derivation of the distribution probabilities



There is another way to extract the probabilities from the probability
generating function (1-ϕ + ϕ F_K(t))^M. Note that this is a
finite degree polynomial in t. The coefficient of t^n gives the
probability F[n] of exactly n amplicons. Instead of generating the
full distribution using the -based algorithms above, we can
extract this single value by using Cauchy's residue theorem with a
contour in the complex plane that contains the origin,

    F[n] = 1/2π i∮(1-ϕ + ϕ F_K(t))^M/t^n+1 dt.


If we take the contour to be the unit circle in the complex plane
centred at the origin, we may make a change of variable,
t = exp(- 2π i θ), so that θ ranges over [0,1], and
the contour integral becomes

    F[n] = ∫_0^1  (1-ϕ + ϕ F_K(exp( -2π i θ)))^M  exp( 2π n θ)dθ

We may evaluate this numerically by splitting the range up into a
large number N of equal sized intervals (with size dθ = 1/N)
and using the trapezoidal approximation. Note that the beginning and
the end point values are identical because we are evaluating the
closed circular contour, so that as the approximation we have


    F[n] ≈1/N∑_j=0^N-1 
    (1-ϕ + ϕ F_K(exp( -2π i j/N))^Mexp( 2π nj/N)

 
 This is fully equivalent to using the above if the same
   value of N is used, and hence is exact if N is
   sufficiently large.
 
 The advantage of this formulation is that if only a single value from
 the distribution is required, it may be found without the large
 overhead in computer memory that using the incurs in storing the
 arrays F[] and F[].

 There is another use for this approach. When evaluating a likelihood
 for a peak height of an allele, it may be that no peak is observed at
 the allele position in the , or one is observed but is not above
 the analytic threshold. In such cases we need to find the cumulative
 probability to the threshold.  Suppose that the threshold corresponds
 to n amplicons. Then the cumulative probability is


    ∑_m=0^n F[n]    = ∑_m=0^n1/2π i∮(1-ϕ + ϕ F_K(t))^M/t^m+1 dt
       = 1/2π i∮  (1-ϕ + ϕ F_K(t))^M
                          (
                          ∑_m=0^n 1/t^m+1) dt
       = 1/2π i∮  (1-ϕ + ϕ F_K(t))^M 1 - t^-n-1/t-1 dt


Again, the trapezoidal rule may be used to evaluate this with line
integral on the unit circle in the complex plane.



 §.§ Simulating the process with stutters

 
The paper of <cit.> also included a model for
stutters. There is a slight error in their paper that was pointed out
in <cit.>. Their corrected model is as follows. During an
amplification an amplicon may make a copy of itself, with probability
p_t, or make a stutter copy with probability p_s. Thus neither a
stutter nor an exact copy is made with probability 1-p_t-p_s.  It is
assumed that when a stutter amplicon is made it will make a copy of
itself in each subsequent cycle with probability p_t.  Let n_k
denote the number of amplicons and m_k the number of stutter
amplicons after k cycles. Assuming N_0 starting amplicons sampled
binomially with probability ϕ for amplification, a simulation for
the number of amplicons and stutters may be expressed as in
alg-gillstutt. In this, the draw of J simulates the total
number of new amplicons and stutters produced in the k-th cycle by
the existing amplicons. The J_t draw then samples from these new
products to simulate how many are copies of the target amplicon, with
the remainder being stutters.  The J_s draw simulates how many new
stutter amplicons are produced by currently existing stutter products.


[Simple amplicon simulation with stutters]
  
  
  
  * Sample n_0 from (N_0,ϕ), and set m_0 = 0.
  
  * for k in 1 to K do
    
    
  * Sample J from (n_k-1,p_s+p_t)
    
  * Sample J_t from (J,p_t/(p_s+p_t))
    
  * Sample J_s from (m_k-1,p_s+p_t)
    
  * Set n_k = n_k-1+ J_t
    
  * Set m_k = m_k-1 + (J-J_t) + J_s
    
  
  

In <cit.> a value of p_s = 0.002 was estimated from
experimental data. However, after correcting their error (which was to
use (n_k-1,p_t) for the number of copies of amplicons
generated, and (n_k-1,p_s) for the number of stutter
products generated from amplicons, hence allowing the possibility for
an amplicon to generate both a copy of itself and also a stutter in a
single cycle) it appears that a value of around 0.004 would be more
appropriate.

 
fig:gillscatter shows a scatterplot of 10000 simulated
(n_K, m_K) values for K cycles starting from a single amplicon,
n_0=1 (so no pre-sampling) using p_t = 0.8 and p_s = 0.004. The
plot indicates some correlation between the two values — for the
data in the plot the correlation coefficient is around 0.22, a
theoretical value will be given later.  <cit.>
also found such correlations in their analysis.




From the simulated data we may also estimate the marginal distribution
for stutter amplicons, a density estimate is shown in
fig:ampliconStutterMarg.
 



We may also find the fraction of the number of amplicons generated
that are stutters, that is the ratio of the number of stutter
amplicons to the total number of amplicon and stutter amplicons, and
generate a kernel density estimate. (The values of this ratio are
therefore between 0 and 1, unlike the stutter proportion which is the
ratio of stutter amplicons to amplicons, which can grow very large.) A
kernel density estimate is shown in fig:gillsratio.






 §.§ Generating the full joint distribution


We now show how to find the full joint- distribution of amplicon
copies and stutter amplicons, at least theoretically.  To do this we
combine a 2-dimensional analysis with bivariate probability
generating functions. A bivariate  is appropriate here, as we
have both main amplicons and stutter amplicons, that is in the
amplification cycles we have a multi-type branching process of the
sort that was formulated by <cit.> in terms of vectorial
generating functions.

We shall refer to the main, or initial, type of amplicon as a
target amplicon, and use the symbol t in the bivariate to represent their number. We shall use the symbol s to represent
the stutter amplicon number. As for the simple (no-stutter) model, we
may express the by functional iteration. However, because we have
two types of amplicons, we need to use coupled equations to express
the amplifications of single target and single stutter amplicons:


    t    →  (1-p)t + p(1-ξ)t^2 + pξ t s
    
      s    → (1-p)s + ps^2
 

Here we let p denote the probability that an amplicon (of either
type) makes a product in a cycle; ξ is the conditional probability
that a target amplicon produces a stutter given either a target or
stutter is produced.  In terms of the previous notation we have
p_s = pξ and p_t = p(1-ξ).  Let F_0(s,t) = t and
G_0(s) = s. Let F_n(t,s) denote the bivariate of the joint
distribution for the number of target and stutter amplicons after n
amplification cycles, arising from a single initial target amplicon,
and let G_n(s) denote the number of stutter amplicons that (would)
arise from an initial single stutter amplicon.

Then for n>0 we have the iteration scheme:


    F_n(t,s)    = (1-p)F_n-1(t,s) + p(1-ξ)F_n-1^2(t,s) + pξ F_n-1(t,s)G_n-1(s) 
    
      G_n(s)    = (1-p)G_n-1(s) + p G_n-1^2(s).


To see this is the case we argue as follows (the reader may find it
helpful to draw a probability tree). Initially we have one target
allele with F_0(s,t) = t. In the first cycle there are three
possible outcomes.

In the first outcome, with probability (1-p) we still have just the
one target amplicon. Hence with a further n-1 cycles the bivariate
for the number target and stutter amplicons conditional on this
outcome is F_n-1(t,s).

The second possible outcome from the first cycle is that the original
target makes a target copy.  This happens with probability p(1-ξ),
and the conditional on this outcome is t^2 representing the two
targets. Now in the subsequent (n-1) cycles each of these targets
will independently give rise to set of target and stutter amplicons.
The for the descendants from each target will each be
F_n-1(t,s). Because the amplifications arising from each of the
two targets are independent, the representing the total number of
target and stutter amplicons from these two targets will be the
product of each of their s, that is F_n-1^2(t,s).

The third outcome is that the initial target produces a stutter. This
happens with probability pξ, and we have the joint conditional on this outcome is st. In the subsequent n-1
amplifications the number of target and stutter amplicons generated by
the target will have F_n-1(t,s). Independently the stutter
amplicon will generate further stutter amplicons with G_n-1(s). The bivariate for the total number of target and
stutter amplicons generated from the single target and amplicon is,
because of the independence of amplification, the product of their
s, that is F_n-1(t,s)G_n-1(s).

Adding these three possibilities together with their probability
weights yields the coupled  of (<ref>): note that the
recurrence relation for G_n(s) is that presented in
thm:pgfs for the non-stutter amplification process.


Now F_n(t,s) given above is for a single starting amplicon. If we
start with M amplicons that are pre-sampled binomially with
probability ϕ, then the final distribution is given by


    ( 1 - ϕ + ϕ F_n(t,s))^M




 §.§ implementation of target and stutter distribution


We now show how the recurrence relations (<ref>) and
(<ref>) may be evaluated numerically using , at least in
principle. (In practice the memory requirements to carry this out will
be excessive for the number of cycles used in forensic analyses.)
The extension is to use a 2-dimensional
. alg-ampliconjointfft gives the details (a Python
implementation is given in py:alg-ampliconjointfft).


[Joint distribution for target and stutter amplicons
    ]
  
  
  
  * Set N = M 2^K
  
  * Initialize F[,] to be a two dimensional N× N array,
    (with lowest index [0,0]) initialized such that all entries are
    zero except F[1,0] = 1.
  
  * Initialize G[] to be an N dimensional array such that all
    entries are zero except G[1] = 1.
  
  * Set F[,] equal to the 2-dimensional of F[,] .
  
  * Set G[] denote the one-dimensional of G[] .

  
  * For K times do:
    
    
  * for each element F[f,g] of F[f,g]
      
      
  * update
        F[f,g]:= (1-p) F[f,g] + p(1-ξ)  F[f,g]^2
        + pξ F[f,g].
      
    
  * for each element G[g] of G[] update
      G[g]:= (1-p) G[g] + p  G[g]^2.
    
  
  * for each element (f,g) of F[] update
    F[f,g]:= (1 - ϕ + ϕ F[f,g])^M
  
  * Set F[,] equal to the inverse of F[].
  
  


Although the use of the alg-ampliconjointfft is not practical
for large K and or M, it can be used for small K and M values.
fig:bivaramplicon shows a contour plot of the joint
distribution for a single amplicon with subject to K=13
amplification cycles, with p=0.85 and ξ = 0.005 on each cycle,
and fig:bivarampsurface a surface plot.  From the figures we
can see quite clearly the correlation between the number of target and
stutter amplicons.  Multi-modality of the distribution is also clearly
evident. The reader may care to compare these figures to
fig:gillscatter.








 §.§ Stutter marginal distribution


From the joint F(t,s) of target and stutter amplicons, we may
obtain the marginal distribution of the target amplicon simply by
substituting s=1. Similarly we may obtain the of the marginal
distribution of the number of stutter amplicons by substituting
t=1. It is the latter we are interested in here. Recall that the
iterative equations finding the joint are


    F_n(t,s)    = (1-p)F_n-1(t,s) + p(1-ξ)F_n-1^2(t,s) + pξ F_n-1(t,s)G_n-1(s) 
    
      G_n(s)    = (1-p)G_n-1(s) + p G_n-1^2(s).

with initial values F_0(t,s) = t and G_0(s) = s. Rather than
iterate to find F_n(t,s) and then substituting t=1 to find the
stutter marginal, we may instead first make the substitution and then
iterate: the equations then simplify to


    F_n(s)    = (1-p)F_n-1(s) + p(1-ξ)F_n-1^2(s) + pξ F_n-1(s)G_n-1(s) 
    
      G_n(s)    = (1-p)G_n-1(s) + p G_n-1^2(s).

with initial values F_0(s) = F_0(1,s) = 1 and G_0(s) = s.

These can be evaluated numerically using a pair of 1-dimensional
, the algorithm is almost identical to
alg-ampliconjointfft, (a Python implementation is given in
py:alg-ampliconstuttermargfft).

[Joint distribution for target and stutter amplicons
    ]
  
  
  
  * Set N = M 2^K
  
  * Initialize F[] to be an N dimensional array such that all
    entries are zero except F[0] = 1.
  
  * Initialize G[] to be an N dimensional array such that all
    entries are zero except G[1] = 1.
  
  * Set F[,] equal to the of F[,] .
  
  * Set G[] equal to the of G[] .

  
  * for each index g of F[]
    
    
  * For K times do:
      
      
  * update
        F[g]:= (1-p) F[g] + p(1-ξ)  F[g]^2 +
        pξ F[g] G[g]
      
  * update
        G[g]:= (1-p) G[g] + p  G[g]^2.
      
      update F[g]:= (1 - ϕ + ϕ F[g])^M
    
  
  * Set F[] equal to the inverse of F[].
  
  

fig:ampliconstuttermarg shows the marginal distribution of
the number of stutter amplicons arising from a single initial target
amplicon, and amplified for K=22 cycles with p=0.8 and
ξ = 0.004 on each cycle. The long right-hand tail of the
distribution is evident.



fig:stutterM83K18 shows the marginal distribution of the
number of stutter amplicons arising from an initial set of 83 target
amplicons, sampled binomially with selection probability ϕ = 2/11
, and amplified for K=18 cycles with p=0.8 and ξ = 0.004 on
each cycle.









 §.§ Moment analysis


The previous sections have shown how the full joint distribution of
target and stutter amplicons may be found numerically by evaluating
the the . Up to now, probabilistic genotyping software
packages, lacking this evaluative ability, assume some simple
distributional assumption— the most common distributions used are
the normal, lognormal of gamma. However, it is straightforward to
derive, from the , expressions for the means and variances of the
full distribution. These moments may be used for finding a `best
fitting' simple distribution based on moment matching.

Let T denote the number of target amplicons, and S the number of
stutter amplicons, arising from the amplification of a single target
amplicon. We have that

    T n    = ∂ F_n(t,s)/∂ t|_t=1,s=1
     T(T-1) n    = ∂^2 F_n(t,s)/∂ t^2|_t=1,s=1
     S n    = ∂ F_n(t,s)/∂ s|_t=1,s=1
     S(S-1) n    = ∂^2 F_n(t,s)/∂ s^2|_t=1,s=1
     TS n    = ∂^2 F_n(t,s)/∂ s∂ t|_t=1,s=1

from which the variances T and T and the correlation
(N,M) may be found.  (Derivations of algebraic solutions are
given in app:ampmoments using different notation.) We have
(variance formulae are given in the appendix)


    T n      = (1 + p(1-ξ))^n
     S n    =  (1+p)^n - ( 1 + p(1-ξ))^n
     TS n    =(1 + p*(1-x))^n-1( [1 - (1+p)^n]x + 2*[(1+p)^n - (1+p*(1-x))^n])
      
     (T,S n)    =  ( TS n) - ( T n )( S n )
       = (1 + p(1-ξ))^n-1[ (1-p(1-ξ))[ (1+p)^n - (1 + p(1-ξ))^n] - ξ((1+p)^n-1)]




As an alternative to using the analytic formula, instead one could fix
numerical values for p, ξ and K and solve the recurrence
relations numerically and very simply.  It is interesting to examine
how the correlation varies with p and ξ and the number of
cycles.  In fig:corr28 we show the dependence for 28 and 34
cycles, in the form of contour plots.

For forensic applications we expect that p will be in the range
0.75 - 0.95, with ξ in the range 0.004 - 0.01. Thus we see
that there can be quite high correlation between the number of
amplicons (which is proportional to the peak heights in the ) of
stutter and target alleles. This is apparently at variance with the
experimental results of <cit.>— we shall
return to this point later.





The moment formulae above assume a single initial amplicon, but they are readily
extended to take account of binomial pre-sampling. There are two ways
to proceed. One way is to write down the to include the sampling;
this is given by


    (1-ϕ + ϕ  F_n(t,s))^M.

This can then be differentiated to find the moments.

Another, simpler way, is to use conditional expectation:

    E Y    = E[E Y X]
    
      V Y    = E[V Y X] + V[E Y X]

for the means and variances of T and S, where in our case
X ∼(M,ϕ) having mean Mϕ and variance
Mϕ(1-ϕ).  Let T denote the total number of amplicons
arising from pre-sampling M target amplicons and then amplifying the
sample, and let S denote the number of stutter amplicon
products.

With this we obtain


    T   = Mϕ [T n] 
     T   = Mϕ ([T n] + (1-ϕ) E[T n]^2)
     S   = Mϕ [S n] 
     S   = Mϕ ([S n] + (1-ϕ) E[S n]^2)


and for the covariance we have

    ( T, S) = Mϕ( (T,S n) +(1-ϕ) [T
      n] [S n])


from which we deduce that the correlation between T and S does not
depend on the initial number of target amplicons, but does depend on
the sampling probability ϕ:

    ( T, S) (T,S n) + (1-ϕ)( T n)( S n))/√(([T n] + (1-ϕ) ( [T n])^2)√(([S n]+ (1-ϕ) ( [S n])^2)



In many forensic laboratory experiments DNA samples are prepared by
high dilution of large template DNA.  For such scenarios a Poisson
distribution would appear to be more appropriate.  Let us assume that
the total number of target amplicons selected for has a Poisson
distribution with rate λ with λ >0.  If
X∼(λ), then the of X has the form


    exp(λ(t-1))


To obtain the for the total number of amplicons we simply
substitute F_n(t,s) for t in this Poisson 
    exp(λ(F_n(t,s)-1))

and take derivatives as appropriate, Alternatively, we may use
conditional expectation noting the mean and variance of a
Poisson(λ) random variable is λ. Either way we obtain


    ( T, S) =  [TS n]/√(([T n]^2)([S n]^2))


Surprisingly, the correlation ( T, S) does not
depend on λ (see app:ampmoments).  Note also that
(<ref>) reduces to (<ref>) with ϕ= 0, which
corresponds to the double limit M→∞, ϕ→ 0 with
Mϕ = λ fixed.

The correlations can get very high, for example with p=0.8,
ξ = 0.005 and k=28 cycles, the correlation is approximately
0.74.






§ AMPLIFYING GENOMIC STRANDS: A GENOMIC MODEL

 
The previous section examined the simple amplicon model of
<cit.>, and showed how by using and the is
it possible to extract, in a very simple manner, the full probability
distributions for the number of target and stutter amplicons. Two
simplifications in the model are that it starts from amplicons, and
not genomic strands, and it does not take into account the dye-tagging
of the amplicons. We shall now remove these limitations; we shall also
add in the artefacts of forward and double stutter, drop-in and
baseline noise, and show how they may all be expressed in a unified
manner using s, and evaluated using s. We being by
considering the genomic strand model without stutters.




 §.§ The basic model described


The basic difference between the model developed here and the
simplified model of <cit.> is that the latter assumed
that we start the branching process with amplicons; in reality we
start the branching process with samples of genomic strands. As
pointed out by <cit.>, this means that
amplicons are not formed until at least the end of the second thermal
amplification cycle as we show in sec:mathformnostutt.



Consider starting from a single genomic strand, and let us assume that
amplification is 100% efficient in all stages. We are interested in a
particular locus.



  * In the first thermal cycle, the genomic strand is melted, and
  primers attach to the flanking regions, one on each of the two
  strands of the melted genome. The TAQ polymerase lays down dNTRs
  along each strand from the primer to past the complimentary flanking
  region. Call these generated sequences half-genomes.  Then
  on cooling we have two separate strands, each consisting of one of
  the complementary strands making up the genomic strand bound to a
  half-genome.

  * In the second thermal cycle these two hybrid genomes melt, and
  we have the two parts of the genome and also the two half-genome
  strands in solution. In this cycle the genomic strands behave as in
  the first thermal cycle, but the half-genome strands, which were not
  there then, attach primers to their complimentary flanking regions,
  and then the TAQ lays down dNTRs along it all the way to the (far)
  end of original flanking region, where the process stops because 
  that is the end of the molecule. Each half-genome strand now has
   a complementary amplicon strand attached to it.

  * In the third thermal cycle the two complementary amplicon
  strands separate from the half-genomes and make complementary copies
  of themselves, and so the exponential growth in the number of
  amplicons with further thermal cycles begins.


The above description is not quite complete, again as described by
<cit.>, in that in addition to making the
amplicons, the amplicons have attached to them a fluorescent dye at
one of the flanking regions, say the 5' region (so that it can be seen
in the capillary electrophoresis (CE) equipment). When measuring the
amplicons in the CE equipment, the amplicons are heated to separate
them out into individual complementary strands, only one each of each
complement has a dye attached.

<cit.> included initial genomic strands and
amplicons in their analysis of , but did not include the
intermediate half-strands. They also did not consider direct double
stutter or forward stutter, or the tagging of amplicons with
fluorescent dyes.



 §.§ Initial mathematical formulation : no stutters



We may denote the original double-helix genomic strand by the pair g
and g_d The half-genomic strands may be denoted by h and h_d,
and the amplicons by a and a_d.  The d subscript denotes a
florescence dye attached. (The original genomic DNA does not have a
dye attached to it, but we use g_d to avoid introducing a special
notation; the half-genomes and amplicons do have such dyes attached.)

The branching process is then summarized by the following processes on
the various components:

    g    → g, h_d
    
      g_d    → g_d, h
    
      h    → h , a_d
    
      h_d    → h_d, a
    
      a    → a, a_d
    
      a_d   → a_d, a

and we are interested in the final number of tagged amplicons, a_d
after the n cycles.

It is also simpler to break the process up into two independent
branching process, one that result in products from the g strand,
and one from the g strand:

    g   → g, h_d
    
      h_d    → h_d, a
    
      a    → a, a_d
    
      a_d    → a_d, a

and from the g_d strand:

    g_d   → g_d, h 
    
      h    → h, a_d
    
      a_d    → a_d, a
    
      a    → a, a_d



The following two tables show how the numbers of each type increase in
each cycle for these two processes treated independently, assuming
100% efficiency in each cycle. The second table follows the pattern
of the first.



  
    cycle    g    h_d    a    a_d 
 0    1    0    0    0 

    1    1    1    0    0 

    2    1    2    1    0 

    3    1    3    3    1 

    4    1    4    7    4 

    5    1    5    15    11 

    6    1    6    31    26 

    7    1    7    63    57 

    8    1    8    127    120 

    9    1    9    255    247 

    10    1    10    511    502 
 






  
    cycle    g_d    h    a_d    a 
 0    1    0    0    0 

    1    1    1    0    0 

    2    1    2    1    0 

    3    1    3    3    1 

    4    1    4    7    4 

    5    1    5    15    11 

    6    1    6    31    26 

    7    1    7    63    57 

    8    1    8    127    120 

    9    1    9    255    247 

    10    1    10    511    502 
 



Of interest is the number of tagged amplicons, which is obtained from
adding up the row entries for a_d in each of the two tables above:




  
    cycle n    a_d 
 0    0 

    1    0 

    2    1 

    3    4 

    4    11 

    5    26 

    6    57 

    7    120 

    8    247 

    9    502 

    10    1013 

  



It is readily verified numerically that these totals are the numbers
2^n-n-1 for each number of cycles n.  (They form the sequence
A000225, called the Eulerian numbers, in the
https://oeis.org/The On-Line Encyclopedia of Integer
  Sequences.)  We can verify this algebraically by induction as
follows. Let us consider the g sequence. It clear that after k
cycles the number of h_d is equal to k, as the g has had k
opportunities to make and h_d.  We see that the values in the third
column are (for k>0) the integers 2^k-1-1, and those in the
fourth column are the Eulerian numbers k reduced by 1, that is
2^k-1-k. We can see this is true for all values given in the
table, we take this as a starting position on an inductive proof for
general k>10.
 
Thus at the k-th cycle, assume that there are n_h_d:k = k copies
of h_d, n_a:k = 2^k-1 copies of a and n_a_d:k = 2^k-1-k
copies of a_d. Then in the next cycle:



  * The number of h_d half-strands increases by 1, as the g
  makes a new copy and the previous copies are not destroyed,
  hence 
    n_h_d:k+1 = n_h_d:k +1 = k+1.


  * The number of a amplicons increases by the number of h_d and
  a_d products that were present on the previous cycle:

    n_a:k+1 = n_a:k + n_h_d:k + n_a_d:k = 2^k-1-1 + k + 2^k-1-k = 2^k-1


  * The number of a_d amplicons increases by the number of a
  amplicons that were present on the previous cycle:

    n_a_d:k+1 = n_a_d:k + n_a_d:k = 2^k-1-k +  2^k-1-1 = 2^k-(k+1)



Hence the inductive hypothesis is proved. Similar calculations go
through for the second table. From this we see that after k+1 cycles
the number of tagged amplicons is


    2^k-1 +  2^k-(k+1) = 2^k+1 - (k+1) -1,


which are the Eulerian numbers.  Note also that the simple amplicon
model would give 2^k+1.  The above formulae correct the values in
Table 4.1 of <cit.>, who ignored the
persistent presence of the original genomic strands and the half
strands in the remaining cycles, assuming their contribution in
subsequent cycles would be negligible and concluding that after n>2
cycles there would be approximately 2^n-2 amplicons.




 §.§ The basic model: PGF formulation



The basic model can be formulated as a multivariate , derived
using vectorial generating functions.  For this we can consider the
s of the two types of initial genomic strands separately.

We start with the g strand, which has the set of amplification
sequences:

    g   → g, h_d
    
      h_d    → h_d, a
    
      a    → a, a_d
    
      a_d    → a_d, a

We introduce symbols t_g, t_h_d t_a and t_a_d to be used
in the multivariate for the number each type of strand.  For a
single strand of each type, each amplifies in a single cycle according
to the s

    t_g    → (1-p_g)t_g + p_gt_gt_h_d
    
      t_h_d   → (1-p_h_d) t_h_d+ p_h_dt_h_dt_a
    
      t_a   →(1-p_a)  t_a+ p_at_at_a_d
    
      t_a_d   →(1-p_a_d) t_a_d + p_a_dt_a_dt_a

where p_g, p_h_d, p_a, p_a_d are branching process probabilities
for each of the types of strands.

Let G_n(t_g, t_h_d, t_a, t_a_d) denote the joint arising
from amplifying a single g strand for n cycles.

Let H_d;n(t_h_d, t_a, t_a_d) denote the joint arising
from amplifying a single h_d strand for n cycles.

Let A_n( t_a, t_a_d)) denote the joint arising from
amplifying a single a strand for n cycles.

Let A_d;n(t_a, t_a_d) denote the joint arising from
amplifying a single a_d strand for n cycles.

Then these s obey the following recurrence relations, which
follow the pattern of the single strand formulae above:

    G_n+1 (t_g, t_h_d, t_a, t_a_d)    = (1-p_g) G_n(t_g, t_h_d, t_a, t_a_d)+
                                            p_gG_n(t_g, t_h_d, t_a, t_a_d)H_d;n(t_h_d, t_a, t_a_d) ,
    
      H_d;n+1(t_h_d, t_a, t_a_d)     = (1-p_h_d) H_d;n(t_h_d, t_a, t_a_d)+ p_h_dH_d;n(t_h_d, t_a, t_a_d)A_n( t_a, t_a_d) ,
    
      A_n+1( t_a, t_a_d)     = (1-p_a)A_n ( t_a, t_a_d))+ p_aA_n( t_a, t_a_d)A_d;n( t_a, t_a_d),
    
      A_d;n+1 ( t_a, t_a_d)    = (1-p_a_d)A_d;n( t_a, t_a_d) + p_a_dA_d;n( t_a, t_a_d)A_n( t_a, t_a_d),

with initial conditions

    G_0 (t_g, t_h_d, t_a, t_a_d)   = t_g
    
      H_d,0(t_h_d, t_a, t_a_d)    = t_h_d
    
      A_d( t_a, t_a_d)    = t_a
    
      A_d,0( t_a, t_a_d)    = t_a_d

It is the last two that give rise to the exponential growth of
amplicons in the process.  Similar equations arise when starting
from a genomic strand g_d, specifically:

    G_d;n+1(t_g_d, t_h, t_a, t_a_d)     = (1-p_g_d) G_d;n(t_g_d, t_h, t_a, t_a_d) + 
                                                p_g_dG_d;n(t_g_d, t_h, t_a, t_a_d) H_n(t_h, t_a, t_a_d) ,
    
      H_n+1(t_h,(t_a, t_a_d)      = (1-p_h)H_n(t_h, t_a, t_a_d)  + p_hH_n(t_h, t_a, t_a_d) B_d;n( t_a, t_a_d) ,
    
      B_d;n+1( t_a, t_a_d)      = (1-p_a_d) B_d;n( t_a, t_a_d) + p_a_dB_n( t_a, t_a_d) B_d;n( t_a, t_a_d) ,
    
      B_n+1( t_a, t_a_d)      = (1-p_a) B_n( t_a, t_a_d) + p_aB_d;n( t_a, t_a_d) B_n( t_a, t_a_d) ,

with initial conditions

    G_d;0(t_g_d, t_h, t_a, t_a_d)     = t_g_d
    
      H_0 (t_h, t_a, t_a_d)    = t_h
    
      B_0(t_a, t_a_d)     = t_a
    
      B_d,0 (t_a, t_a_d)    = t_a_d

where we introduce functions B and B_d in place of A and A_d
as these are to be considered as iterating (amplifying) independently.
The justification of the recurrence relations is similar to that given
in sec:fullampstutdist for the joint distribution of target
and stutter amplicons in the amplicon model, and is omitted.


In general, this is a multivariate polynomial in six t-parameters,
as well as depending on six branching probabilities.  However, because
the genomic strand itself does not duplicate, we see that G and
G_d are proportional to t_g and t_g_d respectively, which is a
simplification that effectively reduces the number of t parameters
to four.

We are particularly interested in the final number of dye-tagged
amplicons. This is the sum of their number arising from the g and
g_d strands, and therefore (because of independence of the branching
process amplifications) is the coefficient of t_a_d in the product
of the s of each strand:

    F_n(t_g, t_g_d, t_h,t_h_d, t_a, t_a_d) = G_n(t_g, t_h_d, t_a, t_a_d)G_d;n(t_g_d, t_h, t_a, t_a_d)

in which we may set all of the t's except t_a_d to unity to
obtain the marginal for the tagged amplicons, that is:

    F_n(1,1,1,1, t_a_d) = G_n(1,1,1, t_a_d)G_d;n(1,1, 1,t_a_d)




 §.§ Moment analysis



Of interest also are the various moments, in particular the mean and
variance, which may be found from the marginal by differentiation
with respect to t_a_d and then setting this to 1.  Algebraic
derivations maybe found in app:genomicmoments, where it is
shown, for example, that the mean number of tagged amplicons is given
by

    p_gp_h_d/p_ap_a_d(
      √(p_a/p_a_d)(1+√(p_ap_a_d))^n-(1-√(p_ap_a_d))^n/2-np_a
    ) + p_g_dp_h/p_ap_a_d(
      (1+√(p_ap_a_d))^n+(1-√(p_ap_a_d))^n/2-1 )

If p_a = p_a_d = p, this simplifies to:

    p_gp_h_d/p^2( (1+p)^n-(1-p)^n/2-np ) +
    p_g_dp_h/p^2( (1+p)^n+(1-p)^n/2-1 )

and if in addition p_g = p_g_d = p_h = p_h_d = p this reduces
to

    (1+p)^n - np -1.

When p=1 this yields 2^n-n-1, the Eulerian numbers given earlier.


There are far too many parameters to explore the behaviour of the
variance in detail here, however one special case could be of
interest, that in which all the amplification parameters take the same
value.  We can compare the behaviour for the mean and variance to that
of the amplicon model.





We see the mean values will be close, viz (1+p)^n - np -1 compared
to (1+p)^n. Somewhat surprisingly, the variance of the genomic
amplification model turns out to be almost exactly half that of the
Gill amplification model for sufficient number of cycles (the ratio is
0.49997 with N=20 and p=0.7, and 0.50000 if we increase to n=24
cycles) and appears to asymptote to almost 0.5 as the number of cycles
increases.  Thus, although the exact formula for the variance is quite
lengthy, for n ≥ 24 cycles taking the variance formula from the
simplified amplicon model, and dividing the result by 2, appears to
give an excellent approximation to it.




 §.§ Full distribution from vectorial PGF


The numerical derivation of the full distribution follows a similar
pattern to the derivation for amplicon model with stutter. Using a
6-dimensional one could, in principle, derive the full joint
distribution for the six types of strand. However the number of g
and g_d strands does not change in each cycle, thus a 4-dimensional
suffices (for amplification of a single genome). The number of
each of the h and h_d strands increases by at most 1 in each
cycle, with at most n-1 of each after n cycles. The numbers of a
and a_d strands will have a maximum of 2^n-n-1, hence the size of
the required is approximately n^2 2^2n, too large for
forensic applications with n between 27 and 34 cycles.

For M genomic strands pre-sampled binomially with probability
ϕ, the joint is given by


    (1-ϕ + ϕ F_n(t_g, t_g_d, t_h,t_h_d, t_a, t_a_d) )^M



However, obtaining the marginal distribution of tagged amplicons can
be found by setting t_g =t_g_d =t_h=t_h_d=t_a= 1, and using 1-dimensional
s of size at most 2^n.  The algorithm is given in
alg:gemomemarg, with an R implementation in
sec:r-fullgill09.

[Marginal target amplicon distribution for the genomic
    model]
  
  
  
  * Set N = M 2^K
  
  * Initialize G[] ,G_D[] ,H[] , H_D[] and A[] to be an
    N dimensional arrays a such that all entries are zero except the
    first, eg, G[0] = 1.
  
  * Initialize A_D[] to be an N dimensional array such that
    all entries are zero except the second, AD[1] = 1.
  
  * Find the of all the arrays, for example, G[].
  
  * Initialize F[] to be an N dimensional array.

  
  * for each index g of G[]
    
    
  * For K times do:
      
      
  * update
        G[g]:= (1-p_g) G[g] + p_g G[g]
          H_D[g]
      
  * update
        G_D[g]:= (1-p_g_d) G_D[g] + p_g_d
          G_D[g] H[g]
      
  * update
        H[g]:= (1-p_h) H[g] + p_h H[g]
          A_D[g]
      
  * update
        H_D[g]:= (1-p_h_d) H_D[g] + p_h_d
          H_D[g] A[g]
      
  * set a =  A[g] and a_d =  A_D[g]
      
  * update A[g]:= (1-p_a)a+ p_a aa_d
      
  * update A_D[g]:= (1-p_a_d)a_d+ p_a_daa_d
      
      update
      F[g]:= (1 - ϕ + ϕ G[g] G_D[g])^M
    
  
  * Set F[] equal to the inverse of F[].
  
  


fig:genomeamp shows the tagged amplicon distribution for the
genomic strand model, with all amplification probabilities equal to
p=0.9, for K=14 cycles (M= 1 and ϕ = 1); in
red is shown the corresponding amplicon model
distribution



We see that the genomic plot is somewhat smoother than the amplicon
model, and more peaked as well, confirming visually the lower variance
discussed above. The plot on the right shows the results on lowering
p → 0.8. The additional green curve shows a
normal distribution with mean and variance matching the genomic model
- we see the fit is appears to be quite good. However it is the
vertical distances at each x value that are important for evaluating
likelihoods) and in this regard the ratios of values are quite
different from 1 in places. For example in the range [1000,5000] the
ratio of the normal approximation to full distribution value ranges
from 0.777 to 1.07.  Multiplying many such ratios together for many
observed peaks could lead cumulatively to a gross over-estimation or
under-estimation of the likelihood, if there is a systematic bias.






 §.§ Including single step backward stutter


We now include stutter into the genomic model.  For simplicity we do
not consider either forward stutter or double stutter, nor allow a
stutter to stutter, at this stage.  The iterative equations for
amplifying a g-strand without stutter is:


    t_g    → t_g(1-p_g) + p_gt_gt_h_d
    
      t_h_d   → t_h_d(1-p_h_d) + p_h_dt_h_dt_a
    
      t_a   → t_a(1-p_a) + p_at_at_a_d
    
      t_a_d   → t_a_d(1-p_a_d) + p_a_dt_a_dt_a


To include stutter we allow there to be stutter variants of h_d, a
and a_d strands, which we denote with an extra s suffix. We need
extra t's to corresponding to each of these variants, and also
additional amplification probabilities for the new amplification
possibilities.  We could introduce many stutter probabilities as well,
but will keep for simplicity to a single conditional stutter
probability.  Our non-stutter equations become extended to:


    t_g    → t_g(1-p_g) + p_g(1-ξ)t_gt_h_d + p_gξ t_gt_h_sd
    
      t_h_d   → t_h_d(1-p_h_d) + p_h_d(1-ξ)t_h_dt_a+ p_h_dξ t_h_dt_a_s
    
      t_h_sd   → t_h_sd(1-p_h_sd) + p_h_sd t_h_sd  t_a_s
    
      t_a   → t_a(1-p_a) + p_a(1-ξ)t_at_a_d + p_aξ t_at_a_sd
    
      t_a_d   → t_a_d(1-p_a_d) + p_a_d(1-ξ) t_a_dt_a + p_a_dξ t_a_dt_a_s
    
      t_a_s   → t_a_s(1-p_a_s) + p_a_st_a_st_a_sd
    
      t_a_sd   → t_a_sd(1-p_a_sd) + p_a_sdt_a_sdt_a_s


For the initial g_d strand we have to introduce the additional
stutter variant h_s and symbols t_g_d and t_h_s.  The
g_d equations are thus:

    t_g_d   → t_g_d(1-p_g_d) + p_g_d(1-ξ)t_g_dt_h + p_g_dξ t_g_dt_h_s
    
      t_h   → t_h(1-p_h) + p_h(1-ξ)t_ht_a_d+ p_hξ t_ht_a_sd
    
      t_h_s   → t_h_s(1-p_h_s) + p_h_s t_h_s t_a_s
    
      t_a   → t_a(1-p_a) + p_a(1-ξ)t_at_a_d + p_aξ t_at_a_sd
    
      t_a_d   → t_a_d(1-p_a_d) + p_a_d(1-ξ) t_a_dt_a + p_a_dξ t_a_dt_a_s
    
      t_a_s   → t_a_s(1-p_a_s) + p_a_st_a_st_a_sd
    
      t_a_sd   → t_a_sd(1-p_a_sd) + p_a_sdt_a_sdt_a_s


Hence we need to keep track of 10 t's corresponding to the 2 genomic
strands, the 4 half strands and 4 amplicon strands. For the marginal
stutter distribution we are interested in the tagged stutter
amplicons, t_a_sd.

These single-strand equations can be `lifted' to iterative equations
for the , giving generalizations of alg:gemomemarg that
are omitted, for obtaining the marginal distributions of target and
stutter amplicons.  fig:genampstutt shows the marginal
distribution of stutter amplicons for the genomic model, with the red
curve showing the marginal from the amplicon model.




We see that the genomic stutter distribution is narrower and more
peaked, thus having a lower variance than the amplicon model. For the
parameters in the plot, the means and variances are shown in
tab:gneampmv. The means are approximately equal, whereas the
variances are about a factor of 2 different. The derivation of moments
is deferred to app:gmom2. fig:genamp2 shows an
example of a marginal stutter distribution with binomial pre-sampling
of strands.








 §.§ Extension to forward stutter and double stutter


The extension to include forward stutter and double back-stutter is
straightforward, and follows the path of stutter.  We need two extra
h variables for the forward stutter and double stutter products, and
4 extra a variables, two for forward stutter and two for the double
stutter. Thus we have 16 variables altogether.

Let the subscript f denotes forward stutter, and r double reverse
stutter. Let ξ_r denote the conditional probability of double
stutter in a cycle, and ξ_f that for forward stutter.  Then the
individual branching equations are given by:

    t_g    → t_g(1-p_g) + p_g(1- ξ_r -ξ_s - ξ_f )t_gt_h_d+ 
            p_gξ_r t_gt_h_rd + p_gξ_s t_gt_h_sd+
            p_gξ_f t_gt_h_fd
    
      t_g_d   → t_g_d(1-p_g_d) + p_g_d(1- ξ_r -ξ_s -ξ_f)t_g_dt_h +
                p_g_dξ_r t_g_dt_h_r+  p_g_dξ_s t_g_dt_h_s+ p_g_dξ_f t_g_dt_h_f
    
      t_h   → t_h(1-p_h) + p_h(1- ξ_r -ξ_s-ξ_f)t_ht_a_d+
              p_hξ_r t_ht_a_rd+ p_hξ_s t_ht_a_sd+ p_hξ_f t_ht_a_fd
    
      t_h_r   → t_h_r(1-p_h_r) + p_h_r(1-ξ_f)t_h_rt_a_rd
                + p_h_rξ_f t_h_rt_a_sd
    
    
    
      t_h_s   → t_h_s(1-p_h_s) + p_h_s (1-ξ_s-ξ_f)t_h_s t_a_sd
                  +  p_h_sξ_s t_h_s t_a_rd
                  +  p_h_sξ_f t_h_s t_a_d
    
      t_h_f   → t_h_f(1-p_h_f) + p_h_f (1-ξ_s)t_h_f t_a_fd+ p_h_fξ_s t_h_f t_a_d
    
      t_h_d   → t_h_d(1-p_h_d) + p_h_d(1-ξ_r -ξ_s-ξ_f)t_h_dt_a+
                p_h_dξ_r t_h_dt_a_r+ p_h_dξ_s t_h_dt_a_s+ p_h_dξ_f t_h_dt_a_f
    
      t_h_rd   → t_h_rd(1-p_h_rd) + p_h_rd(1-ξ_f)t_h_rdt_a_r
                   + p_h_rdξ_f t_h_rdt_a_s
    
      t_h_sd   → t_h_sd(1-p_h_sd) + p_h_sd(1-ξ_f) t_h_sd  t_a_s + p_h_sdξ_f t_h_sd  t_a
    
      t_h_fd   → t_h_fd(1-p_h_fd) + p_h_fd (1-ξ_s)t_h_fd t_a_f+ p_h_fdξ_s t_h_fd t_a
    
    
      t_a   → t_a(1-p_a) + p_a(1-ξ_s - ξ_f)t_at_a_d + p_aξ_s t_at_a_sd + p_aξ_f t_at_a_fd
    
      t_a_d   → t_a_d(1-p_a_d) + p_a_d(1-ξ_s-ξ_f) t_a_dt_a + p_a_dξ t_a_dt_a_s+ p_a_dξ_f t_a_dt_a_f
    
      t_a_s   → t_a_s(1-p_a_s) + p_a_s(1-ξ_s-ξ_f)t_a_st_a_sd
                + p_a_sξ_s t_a_st_a_rd
                + p_a_sξ_f t_a_st_a_d
    
      t_a_sd   → t_a_sd(1-p_a_sd) + p_a_sd(1-ξ_s-ξ_f)t_a_sdt_a_s
                   + p_a_sdξ_s t_a_sdt_a_r
                   + p_a_sdξ_f t_a_sdt_a
    
      t_a_r   → t_a_r(1-p_a_r) + p_a_r(1-ξ_f)t_a_rt_a_rd
                + p_a_rξ_r t_a_rt_a_sd
    
      t_a_rd   → t_a_rd(1-p_a_rd) + p_a_rd(1-ξ_f)t_a_rdt_a_r
                   + p_a_rdξ_f t_a_rdt_a_s
    
      t_a_f   → t_a_f(1-p_a_f) + p_a_f(1-ξ_s)t_a_ft_a_fd + p_a_fξ_s t_a_ft_a_d
    
      t_a_fd   → t_a_fd(1-p_a_fd) + p_a_fd(1-ξ_s)t_a_fdt_a_f
                   + p_a_fdξ_s t_a_fdt_a



The corresponding recurrence relations for the multivariate in
these 18 variables is left to the reader.

With the introduction of forward and double stutter we have the
possibility that over several cycles a stutter product may amplify by
forward stuttering to make target amplicon, or could itself stutter to
make a double stutter. Similarly a forward stutter product from a
target could amplify and stutter thus producing a target, or could
double stutter and product a stutter. Likewise a double stutter
product could amplify by forward stuttering and create a stutter
product. This means that there is a cyclic dependency in the equations
above, so care is required in their computer implementation.

Note that the above equations do not describe the possibility that a
double-stutter strand could itself stutter or double stutter when
amplified, or a forward stutter strand could itself forward
stutter. These more general possibilities are accounted
in part:newframework, in which we give a general modelling
framework, and also include the additional artefacts of drop-in and
background noise, and DNA degradation, that we have not dealt with so
far. But first we look at some other mathematical computational
aspects of the branching process.




§ FURTHER MATHEMATICAL AND COMPUTATIONAL ASPECTS


In this section, further mathematical properties of the branching
process as revealed by the probability generating functions are
explored.  Readers interested in forensic applications may wish to
skip this part and go onto part:newframework.




 §.§ Finding  marginal distributions without a full 

We consider again the simple amplicon model, in which we start with a
single amplicon with recurrence relation

    F_k+1(t) = (1-p)F_k(t)  + pF_k(t) ^2

and initial condition F_0(t) = t.  If we let N=2^K, then we can
find the full distribution using the , or alternatively, for a
particular value n we can find the probability of exactly n
amplicons via Cauchy's residue theorem


    F[n]    = 1/2π i∮ F_K(t)/t^n+1dt  
       = 1/N∑_j=0^N-1F_K(e^-i2π j/N)e^i2π jn/N

To evaluate this sum, we have to evaluate each term
A_K(e^-i2π j/N) for points located on the unit circle in the
complex plane. This can be done the recurrence equation
(<ref>). If one views (<ref>) as a nonlinear
discrete dynamical system, the we see that it has fixed points at
t=0 and t=1. The fixed point at t=0 is stable, whilst the fixed
point at t=1 is unstable: for real values of t∈ [0,1), iterates
of (<ref>) converges to zero, whilst for real t>1
iterates diverge to infinity. It turns out that, apart from the j=0
term which corresponds to the stable point t=1, the iteration scheme
t → (1-p)t + p*t^2 on the unit circle forms a contraction mapping
so that the F_k(e^-i2π j/N) iterates converge to the origin as
k→∞. This is illustrated in fig:cxampiters. (See
also <cit.>.)



Not only do the F_k(e^-i2π j/N) converge to the origin, they do
so rapidly with increasing j. Hence, we may approximate
(<ref>) by a truncated series with a point symmetrically
located around j=0 (ensuring the sum total is real valued), thus:


    F[n] ≈1/N∑_j=-L^LF_K(e^-i2π j/N)e^i2π jn/N


fig:fulltraunc overlays the exact distribution obtained from
a full analysis, and the distribution for the truncated
approximation of (<ref>). The curves are visually
indistinguishable.





Taking account of binomial pre-sampling with truncation requires a
little care. The full series expansion, for M starting amplicons, is


    1/N∑_j=0^N-1(1-ϕ + ϕ F_K(e^-i2π j/N))^M
    e^i2π jn/N


To carry out the truncation, we set F_K(e^-i2π j/N))=0 on the
j terms on the circle further than ± L from the j=0 term, to
obtain

    1/N∑_j=-L^L(1-ϕ + ϕ F_K(e^-i2π j/N))^M
    e^i2π jn/N + 1/N∑_j=-L*^L*(1-ϕ )^M e^i2π
      jn/N

where the second summation is over terms complementary to the first
summation. We now use that e^i2π n/N is an N-th root of unity,
for which ∑_j=0^N-1e^i2π jn/N = 0, to obtain

    ∑_j=-L*^L*(1-ϕ )^M e^i2π jn/N =
    -∑_j=-L^L(1-ϕ)^M e^i2π jn/N

so that

    P[n] ≈1/N∑_j=-L^L(  
        (1-ϕ + ϕ F_K(e^-i2π j/N)^M - (1-ϕ)^M
      )
      e^i2π jn/N



Drop-in can be handled in a similar manner, as can finding stutter
marginal distribution values.  Obtaining good convergence typically
requires more terms than for (<ref>). Details are left to the
interested reader to investigate (or assign to a research student if
available).  Note that if the series is truncated prematurely it could
lead to poor, and even negative, probabilities.


Bivariate distributions may also be found using similar truncated
series. As pointed out in sec-ampliconjointfftPy, to do a
full analysis to obtain the bivariate distribution of a target
and stutter amplicons numbers would require an impossible amount of
computer memory and time. However, by writing the 2-Dfunction out
as a double series and truncating to a smaller set of terms near the
origin, efficient and accurate evaluation of individual bivariate
probabilities may be carried out. Accuracy can be gauged by seeing how
the values change with using an increasing number of terms in the
double series. The details are left to the reader to fill in:
fig:bivargnupt shows two views of the bivariate distribution
for a single starting amplicon, of the number of target and stutter
amplicons, using summation limits of L = ± 10^12 for each
evaluated point (thus around 2^26≈67 million terms
altogether).









PART: 
New Framework in detail


The previous sections have given some indication of the possible
models that can be formulated, and how they may be analyzed using the
mathematical tools of multivariate s combined with s.  In
this part in sec:modgen we present a general framework of the
process leading to an , or set of s in the case of replicate
analyses.  Specialization of this framework to a specific model is
then presented in sec:particularmodel, and its performance is
illustrated using simulated in sec:simdata and publicly
available real data in sec:budata.






§ MODELLING THE GENERATION PROCESS




 §.§ Contributor DNA


We assume that we have a sample of C DNA cells in a volume V. The
DNA is assumed to come from I individuals, with the i-th
individual denoted by K_i, who is a contributor of c_i cells to
the sample. Thus C = ∑_i=1^I c_i.  Let L denote the set of
STR loci under investigation, and let A_l denote the set of alleles
for a given locus l ∈ L.  The total number of genomic strands of
type a_l ∈ A_l within the cells contributed by K_i depends on
the genotype of the individual. Let n_ia_l denote the number of
alleles of type a_l∈ A_l of individual K_i.  The values that
n_ia_l can take depends on the locus l. With a few exceptions,
the possibilities are as follows.



  * If l is autosomal, then n_ia_l∈{ 0,1,2}.

  * if l is Amelogenin, then n_iX∈{ 1,2} and
  n_iY∈{ 0,1}. If K_i is male then n_iX = n_iY = 1;
  otherwise K_i is female with n_iX = 2 and n_iY = 0.

  * If l is a Y-linked locus, then either
  n_ia_l∈{ 0,1,2} or n_ia_l∈{ 0,1}, depending on
  the precise locus.

  * if If l is an X-linked locus, then n_ia_l∈{ 0,1} if
  K_i is male, and n_ia_l∈{ 0,1,2} if K_i is female.


Let g_a_l and g_a_l:d denote the two complimentary strands for
the genome of type a_l ∈ A_L. We use t_g_a_l and
t_g_a_l:d for symbols in their generating functions. It is not
until the process begins that the g_a_l and g_a_l:d
strands become separated. Prior to they are combined, hence the
of the number of both of these genomic strands in the DNA sample
of allele type a_l from K_i is
(t_g_a_lt_g_a_l:d)^n_ia_lc_i, and the multivariate of all alleles in A_l from person K_i is

    ∏_a_l∈ A_l(t_g_a_lt_g_a_l:d)^n_ia_lc_i.

If, for example, l is autosomal and K_i is homozygous with
genotype (a_l,a_l) then the is
(t_g_a_lt_g_a_l:d)^2c_i.

Using enzymes, DNA is extracted from the cells in the volume V. Let
π_e:i,a_l denote the extraction efficiency, that is, the
probability that a genome segment containing the allele a_l ∈ A_l
from K_i is extracted in a state suitable for amplification.  The
efficiency π_e:i,a_l will depend on the process of the DNA
extraction, and also the allele via its length. The extraction
process, as well as removing genomes from cells, breaks them up into
small pieces. The location of the breakage points along a genome can
be considered a Poisson process, so that the distance between breaks
has an exponential distribution. If a break occurs in or between the
flanking regions of a genomic strand pair, then it cannot be
amplified, hence the dependence of π_e:i,a_l on the
allele.[With the Poisson process just described, i
  π_e:i,a_l will decay exponentially at a rate proportional to
  the total base-pair size of the between and including the flanking
  regions. This is one possible mechanism of preferential
    amplification <cit.>.]  The for
extraction of the (single) pair is thus

    E_x( t_g_a_l,t_g_a_l:d) =  1 - π_e:i,a_l + π_e:i,a_l t_g_a_lt_g_a_l:d .


However, prior to extraction the DNA may be degraded by age,
environmental or other factors. Such degradation could also be person
specific (for example, depending on the cell type of the DNA from the
person).

One form of degradation leads to a break in the genome-pair strand in
or between the flanking regions. If we let λ denote the
probability of there not being such a break, then the for
extraction, of amplify-able (single) genome pair, becomes

    E_x( t_g_a_l,t_g_a_l:d) =  1 - π_e:i,a_lλ + π_e:i,a_lλ t_g_a_lt_g_a_l:d.



Alternatively degradation could lead to one of the complimentary
strand pairs having a break and the other not, as illustrated in
fig:degrade2.  If we assume such breaks occur independently
for each if the complimentary strands, then the joint for
extraction of the (single) genome pair strands becomes

    Ex( t_g_a_l,t_g_a_l:d) = 1 - π_e:i,a_l(1 - (1-λ)(1-λ_d))+ π_e:i,a_l (1-λ +λ t_g_a_l)(1-λ_d + λ_d t_g_a_l:d).

where we allow for the possibility of different breakage
probabilities, λ and λ_d, for the two strands (arising
from their different, but complimentary, chemical base composition).



A fraction π_f of the volume V is taken for amplification. Hence
the joint for successful extraction and selection of the genome
strands from a single pair is

    1 -π_f  + π_f E_x( t_g_a_l,t_g_a_l:d)



The total number of such strands from person K_i extracted and
selected for amplification and are amplify-able therefore has the 
  
    (1 -π_f  + π_f E_x( t_g_a_l,t_g_a_l:d)) ^n_ia_lc_i.


where E_x( t_g_a_l,t_g_a_l:d) depends on the breakage model
chosen.



The full multivariate for the contributors' strands from all loci
under consideration ready for amplification is therefore


    ∏_l∈  L∏_a ∈ A_l∏_i ∈ I(1 -π_f  + π_f E_x( t_g_a_l,t_g_a_l:d) )^n_ia_lc_i.





 §.§ Including Drop-in


In addition to the extracted alleles, there is a possibility that
spurious drop-in alleles may get into the minitube in which
the aliquot is ready for the process.  We follow
<cit.> and assume that such drop-in events occur by a
Poisson process. Denote the drop-in rate for allele a_l ∈ A_l by
λ_a_l. Then the for the total number of genome strand
pairs that drop-in is given by


    D_r( t_g_a_l,t_g_a_l:d) =  exp( λ_a_l (  t_g_a_lt_g_a_l:d-1))


Drop-ins occur independently for each locus, and independently of the
genotypes of the I individuals. Hence the for the genomic
strands from drop-in and contributors prior to amplification is

    ∏_l∈  L∏_a ∈ A_l∏_i ∈ I(1 -π_f  + π_f E_x( t_g_a_l,t_g_a_l:d) )^n_ia_lc_i
      D_r( t_g_a_l,t_g_a_l:d)




Within the genomic strand framework presented here, there is an
alternative model for drop-in, in which amplicon-pair strands, rather
than genomic-pair strands, fall into the minitube that the is
carried out in. One reason amplicons could drop-in 
is that there is a build-up over time of
contaminating amplicons from the previous analyses carried out in
 the forensic
laboratory. This
possibility is supported by the observation that there are lower
levels of drop-in in samples amplified after a laboratory has been
deep-cleaned, compared to just before a deep-cleaning has been carried
out, and that as time goes on dropin-rates increase until the
laboratory is cleaned again[Sue Pope, personal
  communication.]

To model this alternative, we use a Poisson drop-in model for the
number of amplicons that drop-in, given by the 
    D_r( t_a_a_l,t_a_a_l:d)= exp( λ_a_l (  t_a_a_lt_a_a_l:d-1)),

instead of the genomic strand drop-in model (<ref>).
Hence the for the genomic strands and amplicons prior to
amplification is now

    ∏_l∈  L∏_a ∈ A_l∏_i ∈ I(1 -π_f  + π_f E_x( t_g_a_l,t_g_a_l:d) )^n_ia_lc_i
      D_r( t_a_a_l,t_a_a_l:d).


One could even consider both Poisson processes occurring
simultaneously so that we allow drop-in of both genomic strands and
amplicon strands, leading to


    ∏_l∈  L∏_a ∈ A_l∏_i ∈ I(1 -π_f  + π_f E_x( t_g_a_l,t_g_a_l:d) )^n_ia_lc_i
      D_r( t_g_a_l,t_g_a_l:d) D_r( t_a_a_l,t_a_a_l:d).

in which we can allow a different drop-in rate λ_a_l for
the two processes.




 §.§ amplification

During the process we may generate strands of type
g_a_l, g_a_ld, h_a_l, h_a_ld, a_a_l and a_a_ld.  For
any of these strands to make a successful copy (or copy of some
stutter type) a primer has to bind at the primer binding site. The
probability of this happening will depend on various factors in how
the is carried out (e.g., the kit, temperature etc.). It could
also depend on the type of strand considered. For example the g
strands could be quite long and coil up to create a barrier preventing
a primer from binding. It also depends on the binding energy of the
primer to the strand, and each end of the flanking region generally
has different base-pair composition, so we should expect the binding
energies to be different. We thus introduce primer binding
probabilities specific to each complementary-strand

    p_l,g, p_l,gd, p_l,h, p_l,hd, p_l,a p_l,ad.

Binding a primer is a pre-requisite for a successful copy. However
during the remainder of the duplication process a copying error could
occur, leading to a stutter variant. The probability that this
happens, and the resulting stutter artefact, will depend on the
initial allele a_l and the resulting allele b_l.  We thus
introduce ξ_b_l a_l to denote the conditional probability
that a strand of allele type a_l creates strand of allele type b_l
given that some product is produced We have that for all a_l∈ A_l,

    ∑_b_l ∈ A_lξ_b_l a_l = 1


Note that we have the conditional probability for copying without
stuttering, given some sort of copy if made, is

    ξ_a_l a_l = 1 - ∑_b_l ∈ A_l: b_l  a_lξ_b_l
      a_l.

We shall assume that these probabilities remain constant throughout
all of the cycles of the process, though this assumption could be
relaxed.


The set of s for single-strand single-cycle duplication is then
given by

  
    g_a_l   : (1 - p_l,g) t_g_a_l +   p_l,g∑_b_l ∈ A_lξ_b_l a_l t_g_a_lt_h_a_ld
    
        g_a_ld   : (1 - p_l,gd) t_g_a_l +  p_l,gd∑_b_l ∈ A_lξ_b_l a_l t_g_a_ldt_h_a_l
    
        h_a_l   : (1 - p_l,h) t_h_a_l +  p_l,h∑_b_l ∈ A_lξ_b_l a_l t_h_a_lt_a_a_ld
    
        h_a_ld   : (1 - p_l,hd) t_h_a_ld +  p_l,hd∑_b_l ∈ A_lξ_b_l a_l t_h_a_ldt_a_a_l
    
        a_a_l   : (1 - p_l,a) t_a_a_l +  p_l,a∑_b_l ∈ A_lξ_b_l a_l t_a_a_lt_a_a_ld
    
        a_a_ld   : (1 - p_l,ad) t_a_a_ld +   p_l,ad∑_b_l ∈ A_lξ_b_l a_l t_a_a_ldt_a_a_l


These single strand s are `lifted' to the full vectorial of
the process, as follows. Let F_n(𝐭_l s ) denote the
multivariate for all the possible types of strands that are
generated from a single strand of types

    s ∈{ g_a_l, g_a_ld, h_a_l, h_a_ld, a_a_l,a_a_ld: a_l ∈ A_l}.
 The 𝐭_l represents all the possible symbols for all the possible alleles and their strand types (there will be 6  A_l such symbols). Then the recurrence relations for the branching process have the form:


  
    F_n+1(𝐭_l g_a_l)   = (1 - p_l,g) F_n(𝐭_l g_a_l ) + 
                                          p_l,g∑_b_l ∈ A_lξ_b_l a_lF_n(𝐭_l g_a_l ) F_n(𝐭_l h_a_ld)
    
        F_n+1(𝐭_l g_a_ld)    =(1 - p_l,gd) F_n(𝐭_l g_a_ld) +  
                                            p_l,gd∑_b_l ∈ A_lξ_b_l a_l F_n(𝐭_l g_a_ld) F_n(𝐭_l  h_a_l)
    
        F_n+1(𝐭_l h_a_l)    = (1 - p_l,h) F_n(𝐭_l h_a_l) +  
                                           p_l,h∑_b_l ∈ A_lξ_b_l a_l F_n(𝐭_l h_a_l) F_n(𝐭_l a_a_ld) 
    
        F_n+1(𝐭_l h_a_ld)    = (1 - p_l,hd) F_n(𝐭_l h_a_ld) +  
                                            p_l,hd∑_b_l ∈ A_lξ_b_l a_lF_n(𝐭_l h_a_ld) F_n(𝐭_l h_a_l) 
    
        F_n+1(𝐭_l a_a_l)    = (1 - p_l,a) F_n(𝐭_l a_a_l) +  
                                           p_l,a∑_b_l ∈ A_lξ_b_l a_l F_n(𝐭_l a_a_l) F_n(𝐭_l a_a_ld) 
    
        F_n+1(𝐭_l a_a_ld)    = (1 - p_l,ad)F_n(𝐭_l a_a_ld) +  
                                            p_l,ad∑_b_l ∈ A_lξ_b_l a_l F_n(𝐭_l a_a_ld) F_n(𝐭_l a_a_l)


with initial conditions

  
    F_0(𝐭_l g_a_l)=    t_g_a_l
    
        F_0(𝐭_l g_a_ld) =    t_g_a_ld
    
        F_0(𝐭_l h_a_l) =    t_h_a_l
    
        F_0(𝐭_l h_a_ld) =    t_h_a_ld
    
        F_0(𝐭_l a_a_l) =     t_a_a_l
    
        F_0(𝐭_l a_a_ld) =    t_a_a_ld



The joint of all types of strand products after n cycles
generated from a single initial pair strand (g_a_l,g_a_ld) is

    F_n(𝐭_l g_a_l ) F_n(𝐭_l g_a_ld )





 §.§ The joint after 


If using the genomic drop-in model (<ref>), one now
substitutes t_g_a_l→ F_n(𝐭_l g_a_l ) and
t_g_a_ld→ F_n(𝐭_l g_a_ld ) into
(<ref>) to obtain the joint of all types of product:


    ∏_l∈  L∏_a ∈ A_l∏_i ∈ I
      (1 -π_f  + π_f E_x(  F_n(𝐭_l g_a_l )  ,F_n(𝐭_l g_a_ld ) ) )^n_ia_lc_i
      D_r( F_n(𝐭_l g_a_l )F_n(𝐭_l g_a_ld ) )


Alternatively, if using the amplicon drop-in model
(<ref>), one substitutes
t_g_a_l→ F_n(𝐭_l g_a_l ) and
t_g_a_ld→ F_n(𝐭_l g_a_ld ) into
(<ref>) together with the additional substitutions of
t_a_a_l→ F_n(𝐭_l a_a_l ) and
t_a_a_ld→ F_n(𝐭_l a_a_ld ) :


    ∏_l∈  L∏_a ∈ A_l∏_i ∈ I
      (1 -π_f  + π_f E_x(  F_n(𝐭_l g_a_l )  ,F_n(𝐭_l g_a_ld ) ) )^n_ia_lc_i
      D_r( F_n(𝐭_l a_a_l )F_n(𝐭_l a_a_ld ) )



Depending on the drop-in model chosen, either (<ref>) or
(<ref>) is the multivariate for all the possible
products. For the capillary electrophoresis we require only the
multivariate for the tagged amplicons. This is obtained by
setting all of the 𝐭_l components to 1 except for the
symbols representing the tagged amplicons, that is, the set
{t_a_a_ld: l∈ L, a_l∈ A_l}.  Let 𝐭_ld denote
𝐭_l with this substitution made. Then (<ref>)
with this substitution becomes


    ∏_l∈  L∏_a ∈ A_l∏_i ∈ I
      (1 -π_f  + π_f E_x(  F_n(𝐭_l g_a_l )  ,F_n(𝐭_ld g_a_ld ) ) )^n_ia_lc_i
      D_r( F_n(𝐭_ld g_a_l )F_n(𝐭_l g_a_ld ) )


with a similar equation for substitution into
(<ref>). Hence, taking into account the two variations of
E_x(·, ·) given earlier, we have four modelling
possibilities.



 §.§ RFU scaling and baseline noise

During the capillary electrophoresis phase, a thin tube is dipped into
the amplified product and a high voltage is applied. This forces a
fraction of the product into the tube where it is carried along by a
voltage differential. At a specific point along the capillary tube a
lasers excite the dyes on the tagged amplicons, and the amount that
fluoresces is recorded as the RFU reading. The RFU value is
proportional to the number of tagged amplicons.

There is, therefore, some scale factor which we denote by ρ, that
relates the final number of tagged amplicons of each type to the RFU
reading. The factor will depend upon the machinery, but may also be
expected to be dye dependent with all the other factors constant. Thus
if we consider a single locus l, a given peak height RFU value
r_a_l>0 will correspond to the range
[ρ(r_a_l-1/2),  ρ(r_a_l+1/2)] of tagged amplicons of that
allele, and conversely.

However this does not take account the baseline noise that is usually
present. If we denote the discrete probability distribution of the
baseline noise distribution by η_l for the range [0,w], say,
then we may form the of the distribution as

    ∑_j=0^w η_l[j]z^j

Given the locus, the noise distribution will be independent of the
allele type (it may depend on the dye lane that the locus is in). Thus
for each a_l∈ A_l we may form a“tagged-amplicon equivalent”
noise distribution given by

    η_l,a_l := ∑_j=0^w η_l[j] t_a_a_ld^jρ

These s may then be multiplied into (<ref>) and the
resulting multivariate used for assessing peak height likelihoods
for observed data.

Alternatively one could derive peak height distributions from
(<ref>) and convolute them with the baseline noise
distributions, and use the convolved distribution to evaluate peak
height likelihoods.



 §.§ Multiple replicates from a sample


Sometimes more than one amplification is carried out on a sample
from which DNA has been extracted, each amplified sub-sample is called
a replicate.  This maybe modelled using a multinomial as
follows. Let there be R replicates, with π_f_r the fraction of
the sample used for replicate r.

If we ignore breakage through degradation, then for a single genomic
pair strand g_a_lg_a_ld in the sample, the multivariate for
the number of strands (0 or 1) of the pair being in each replicate is

    1 - ∑_r=1^Rπ_f_r + ∑_r=1^Rπ_f_r
    t_g_a_l:rt_g_a_ld:r

where we introduce new symbols to represent the specific replicate the
pair get selected for, indicated by the additional `:r' subscripts.
New symbols for half-strand and amplicons are also required to specify
the replicate they belong to.

To take into account the extraction efficiency, we multiply each
π_f_r by π_e:i,a_l:

    (1 - π_e:i,a_l∑_r=1^Rπ_f_r)  + π_e:i,a_l∑_r=1^Rπ_f_r t_g_a_l:rt_d_a_ld:r  .


To take account of degradation breakage, replace
t_g_a_l:rt_h_a_ld:r by E_x( t_g_a_l:r,t_g_a_ld:r)
in (<ref>):


    (1 - π_e:i,a_l∑_r=1^Rπ_f_r)  + π_e:i,a_l∑_r=1^Rπ_f_rE_x( t_g_a_l:r,t_g_a_ld:r)  .

so that for person K_i,

    ((1 - π_e:i,a_l∑_r=1^Rπ_f_r)  + π_e:i,a_l∑_r=1^Rπ_f_rE_x( t_g_a_l:r,t_g_a_ld:r))^n_ia_lc_i  .






 §.§ Untyped contributors


Equation (<ref>) assumes that the genotypes of all
contributors are known. For a contributor U_j whose genotype is not
known, we replace (<ref>), for autosomal loci, by the 
    ( ∑_a_l∈ A_l p_a_l (t_g_a_lt_g_a_l:d)^c_j )^2

where p_a_l is the allele frequency of allele a_l in the locus
l for the population that U_j comes from.  For a set of J
unrelated untyped individuals we have

    ∏_j=1^J ( ∑_a_l∈ A_l p_a_l (t_g_a_lt_g_a_l:d)^c_j )^2.


This represents intact genomic pair strands in the sample prior to
extraction and without degradation or splitting a fraction into
r ≥ 1 replicates. To take all these into consideration we replace
the t_g_a_lt_g_a_l:d product pairs thus:


    ∏_j=1^J ( ∑_a_l∈ A_l p_a_l 
        ((1 - π_e:i,a_l∑_r=1^Rπ_f_r)  + π_e:i,a_l∑_r=1^Rπ_f_rE_x( t_g_a_l:r,t_g_a_ld:r))
        ^c_j )^2.





 §.§ Another variation of the framework


<cit.> proposed a method for improving the
detection of alleles as described in the Method Summary of their
paper:

  DNA template is first divided into two aliquots. One aliquot is used
  as template for a using a primer mix containing all forward
  primers for loci targeted in the PowerPlex ESI 16 kit (Promega),
  while the other aliquot is amplified with all reverse
  primers. Amplification products are then pooled for use as template
  in a standard with the STR kit primer mix. The forward and
  reverse primer reactions result in a linear amplification of the
  target sequences to boost the amount of template available for ,
  thus reducing the stochastic effects commonly seen with low template
  DNA analysis.


The framework presented so far does not apply to their experimental
set-up; however it is readily adapted to it as follows.  Assume that
the DNA template is split into two equal aliquot parts for the first
κ pre-cycles, and then put back together for n amplifications
cycles.  Consider a single genome-pair (g ,g_d) (dropping the a_l
indexing from earlier notation for simplicity) in the aliquot prior to
splitting the sample.  Let S and S_d denote the two split samples.
If the sample is split into two equal volumes, then the (g,g_d) pair
will be in S with probability 0.5, and S_d with probability 0.5.
Suppose that the sample S has primers that bind to the g strand;
then the number of h_d strand generated by the κ pre-cycles
will be binomially distributed as (κ, p_g) if the
stand pair is in S This is because the g strand produces at most
one h_d strand per cycle. Similarly the number of h strands
generated in S_d in the κ pre-cycles will be binomially
distributed as (κ, p_g_d) if the strand pair is in
S_d.

Hence when the two samples are re-combined for the main n cycles that generate the amplicons, the for the starting number
of g, g_d, h and h_d strands will be


    t_g t_g_d(1+t_h)^κ + (1+t_h_d)^κ/2


One then replaces t_g→ F_n(𝐭 t_g), and so on for the
joint after the n cycles of . Other factors such as an
unequal splitting of the DNA template, extraction efficiency,
degradation, multiple replicates, number of cells and genotypes of
contributors, and modification to the drop-in model, can be taken into
account in a straightforward manner. For example, for the genomic
drop-in model we make an adjustment from
  
    D_r( t_g,t_g_d) =  exp( λ (  t_g t_g_d-1))

  to
  
    D_r( t_g,t_g_d) =  exp( λ(  t_gt_g_d(1+t_h)^κ + (1+t_h_d)^κ/2-1)).




  

§ A PARTICULAR MODEL REALISATION

  

  In this section we present a model that specialises the framework
  given earlier, and which introduces approximations that make the
  model computationally tractable. A satisfying feature of this
  approach is that the mathematical nature of the model approximations
  are clearly specified, and they can be judged on their merits. The
  author has implemented the model in a computationally efficient and
  accurate system, and this implementation is used in the performance
  analysis of the model applied to real and simulated data presented
  below.




 §.§ Model assumptions and approximations

In the model presented we shall assume that we have a single DNA
sample from which a single replicate has been produced.  First the
sample is described, then the assumptions regarding the process
are given, and how the tagged amplicon number is interpreted as measurement in terms of units. We shall assume that the kit being
used has autosomal loci, and perhaps also Amelogenin but no other
sex-linked loci are among the loci of the kit. We also assume that all
contributors are unrelated.



  §.§.§ The sample


We shall assume that we have a mini-tube of volume V_s, that
contains the extracted DNA of I contributors, in which the ith
contributor has contributed c_i cells.  The DNA has been extracted
with an efficiency ψ∈ [0,1], and may have degradation
characterised by a parameter δ which has units of inverse
length (measured in base pairs).  The values of the cell amounts c_i
and the degradation parameter are taken as unknown, and to be
estimated from the .  The extraction efficiency is taken as known.

A very small amount of the sample in V_s is taken to quantify the
concentration of DNA in V_s, and based on this a fraction fraction
π_f of the volume V_s is taken is put into a mini-tube of volume
V, to which is added primers . is carried out with this
aliquot.

Note that quantification of DNA is usually carried out using qPCR, so
that the amount of DNA estimated to be in V is the amplify-able
amount of DNA. If the estimated amount of DNA in V is γ, then
an estimate of the amount of DNA in the sample volume V_s is given
by γ_s = γ/(ψπ_f), if we ignore degradation.


We shall assume that degradation has the effect of breaking a pair
strand in two or more pieces, according to a Poisson process.  For a
given locus l ∈ L, the total number of genomic pair strands of an
allele a_l ∈ A_L, having total size in base-pairs between and
including flanking regions denoted by b_l, that are in the aliquot
volume V and are amplify-able is given by the binomial distribution,
 taken to be

    ( ∑ _i=1^I n_ia_lc_i, ψπ_f exp(-δ b_l)).


We shall assume the genomic drop-in model, with a Poisson model for
the number of drop-in genome strands with drop-rates for every allele
assumed known. We ignore stutter products that may be generated during
the process by drop-in strands.




  §.§.§ The 

We shall use the genomic model described earlier. We assume that all
amplification probabilities, conditional and unconditional, are known
for all allele and loci.

We shall assume that in an amplification cycle, forward stutter,
no-stutter, single reverse stutter and double reverse stutter products
may be directly produced from a strand of type of the allele
a_l. Note that we allow triple stutters arising indirectly via two
amplification cycles. Some alleles might not form one or more of these
stutter products, for example because they are at the low or high end
of the allelic range of A_l; additionally Amelogenin is assumed not
to form stutters.

We only allow stuttering for multiples of the base-pair repeat size of
each locus. Thus while we include a stutter of 9.3→ 8.3, for TH01,
we exclude the stutter possibility 9.3→ 9; this is purely for
computational efficiency, to make the likelihood evaluation tractable.



  §.§.§ Converting amplicon numbers to peak height s


At the conclusion of the , the volume V will have a large number
of tagged amplicons for the various alleles in the amplification
kit. A small number of these are drawn up electrostatically into the
capillary-electrophoresis machine. We shall assume that the number of
tagged amplicons of a specific allelic type drawn up is proportional
to their number (the proportionality depends on other factors such as
the voltage applied; see <cit.> pp.144-145),
and that the peak height is proportional to the number drawn up,
and that the proportionality constant may depend upon the particular
dye.  Hence we introduce dye-lane proportionality factors so that the
peak height generated by the tagged amplicons is proportional to
the total number in the volume V.

To this must be added a random value arising from the baseline
noise. It is assumed that the prior baseline noise distribution is
known, and may be dye-lane dependent.



  §.§.§ Likelihood evaluation


In order to evaluate the likelihood function requires evaluation of
the s F_n(𝐭_ld g_a_l ) and
F_n(𝐭_ld g_a_ld ), based on the assumed known
amplification probabilities and binomial sampling rates and genome
counts, for all the alleles.  Let us write this as


    ∏_a_l ∈ A_L F_n(𝐭_ld g_a_l )  = 
      ∏_a_l ∈ A_L∏_b_l ∈ A_L ( 1 -ϕ_b_l + ϕ_b_l G(t_b,a_ld))^k_b_l,


where G(t_b,a_ld) is the of the joint distribution of all
allelic products arising from a single genomic strand of type b_l,
ϕ_b_l is the binomial sampling probability for such alleles,
of which there are k_b_l in the sample.

This is computationally intractable for two reasons. The first is in
evaluating the individual s G(t_b,a_ld). We therefore make
the approximation in which each G(t_a_ld) factorizes as


    G(t_a_ld) =  G_-2(t_a_ld) G_-1(t_a_ld) G_0(t_a_ld) G_+1(t_a_ld) 
      = ∏_j=-2^1G_j(t_a_ld)


in which G_-2(t_a_ld) is the of the number of tagged
amplicons in double-stutter position, G_-1(t_a_ld) in stutter
position, G_0(t_a_ld) copies of the target allele, and
G_-1(t_a_ld) forward stutter.


Even with this approximation, the computations are still intractable,
we therefor make the following further approximation:

    ( 1 -ϕ_b_l + ϕ_b_l∏_j=-2^1G_j(t_a_ld)  )^k_b_l→∏_j=-2^1 (1 -ϕ_b_l + ϕ_b_lG_j(t_a_ld) )^k_b_l.

Under these two approximations, the total likelihood, given the
genotypes of contributors, will factorise into a product of functions,
each of which depend on a single allele type.  For a particular allele
a , there will be up to four factors arising from alleles 2 repeats
higher (which double stutter to make a, one repeat higher (which
stutter to form a) , alleles of the same type a, and alleles one
repeat lower (which forward stutter to produce a product).

This needs to be multiplied by s from each allele drop-in ,
and combined with baseline noise distribution on each allele, to get
the marginal peak-height distributions for each allele.  Hence for
each allele we have a product of univariate s that can be
convoluted using the to find the marginal distribution on each
allele, which can then be used directly to find the likelihoods on
each allele given the peak height.

The peak height likelihood from all peaks is then found from these
individual allele likelihoods by multiplication—given the genotypes
of the contributors—and such products averaged using population
genotype probabilities over the possible genotypes of the untyped
contributors.



§ MODEL PERFORMANCE WITH SIMULATED DATA


One of the useful features of the framework developed in this paper is
that, because the model is an idealisation of the DNA extraction and
amplification process, not only can realistic simulated single source
and mixtures s be produced, the simulated data can be analysed
by the model using the same parameters as used in simulating the data.
The performance of the model can then be compared to the `gold
standard' simulated values. Additionally, one can tweak the parameters
of the model so that they differ from those used to simulate the data,
so that robustness of the model can be gauged.  It is also possible to
use unrealistic model parameters to isolate and analyse specific
effects, for use with software validation and regression tests.  All
these can be used to judge the limitations of the model. Here we
present a set of simulations of increasing complexity, beginning with
unrealistic single contributor sample simulations.  All the
simulations use population Caucasian data from <cit.>.
All allele frequencies in a locus were increased so that the minimum
allele count of any was 5. (The total allele count was also increased
to yield population probabilities that add to 1, so this is not quite
the same as the 5/2N adjustment of <cit.>
which does not adjust the normalisation). An Fst value of 0.02 was
used in all analyses.  Following <cit.> we use a
locus-wide drop-in rate of 0.021, with allele specific rates equal to
this multiplied by the relative allele frequency in the
unadjusted population counts. Contributor genotypes were
taken from those in the dataset of contributors from the
PROVEDIt Initiative <cit.>,
which we shall return to in sec:budata when we examine the
model performance on experimental data.


All simulations are based on the Identifiler™ kit, in
which we assume that all strand amplification probabilities are
p_g = p_g_d = p_h = p_h_d = p_a = p_a_d = 0.85 for all loci.



 §.§ Single contributor simulations




  §.§.§ Simulations with no stutter and no noise 

In this set of simulations we take all the conditional stutter
probabilities ξ to be zero on all alleles and loci.  With this
assumption, we will have
G_-2(t_a_ld) = G_-1(t_a_ld) = G_+1(t_a_ld) = 1,
and (<ref>) is no longer an approximation, with only
the j=0 terms surviving on both sides of the equation to give an
equality. We are thus, in these simulations, testing the performance
of the model under the approximation of (<ref>).

We shall simulate data the genotype data of subject
 taken from the PROVEDit dataset, and assume that
the profile is known when estimating cell counts.  We take as further
idealisations that there is no drop-in (the drop-in rate is zero), and
that there is no baseline noise.  We use a factor of 800,000 to scale
post-tagged-amplicon numbers to  values, and 28 cycles for
the amplification.

We use maximum likelihood to estimate the cell counts of the
contributor, a value known at the time of simulation. We can thus
examine the predictive accuracy and variability of the cell count
estimates.

In the first set of simulations, we shall take the degradation
parameter to be zero, and shall assume this value when estimating the
cell counts.  We compare the predictions made using the model based on
the analysis and normal, lognormal and gamma distribution models
based on matching the means and variances of the marginal
distribution.  We take the binomial sampling probability ϕ to be
the same for all loci and alleles, and see how the models behave as we
vary the number of cells, and the values of ϕ. We include in the
simulations the unrealistic value ϕ=1, so that
(<ref>) is exact, and we can compare the (now exact)
model to the moment approximation models.

The following table shows maximized log-likelihoods and estimated
cells counts from 10 simulations, with 500 cells, ϕ =1, and an
analytic threshold of 1 in all simulations; we see that the
estimates are very close to the true value of 500 cells used to
generate the simulations.  We also see that the log-likelihood maxima
are all of a similar value.



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -171.356    499    -171.341    499    -171.345    499    -171.367    499 

    -168.651    499    -168.658    499    -168.655    499    -168.648    499 

    -170.734    500    -170.757    500    -170.749    500    -170.719    500 

    -169.804    499    -169.802    499    -169.802    499    -169.806    499 

    -171.148    500    -171.176    500    -171.166    500    -171.129    500 

    -173.965    500    -174.017    500    -173.999    500    -173.931    500 

    -173.689    498    -173.770    498    -173.743    498    -173.635    498 

    -170.764    499    -170.665    499    -170.697    499    -170.830    499 

    -172.017    500    -172.036    500    -172.029    500    -172.008    500 

    -172.890    499    -172.851    499    -172.863    499    -172.919    499 
 



Reducing the initial number of cells to 10 , we obtain the following
table from 10 simulations, in which every model estimated the actual
number of cells used. Again the log-likelihood values are very similar
in each simulation across the models.


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -110.371    10    -110.321    10    -110.322    10    -110.468    10 

    -106.731    10    -106.708    10    -106.716    10    -106.746    10 

    -116.895    10    -116.288    10    -116.455    10    -117.502    10 

    -105.7    10    -105.645    10    -105.663    10    -105.767    10 

    -107.805    10    -107.955    10    -107.896    10    -107.753    10 

    -111.267    10    -111.747    10    -111.564    10    -111.012    10 

    -112.118    10    -113.147    10    -112.767    10    -111.622    10 

    -110.533    10    -111.237    10    -110.984    10    -110.134    10 

    -106.879    10    -106.954    10    -106.926    10    -106.842    10 

    -119.358    10    -120.664    10    -120.166    10    -118.811    10 

    



We now reduce the ϕ value to 0.7.  For 500 cells, we obtain the
following estimates:




  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -200.822    495    -200.725    495    -200.749    495    -200.853    495 

    -193.161    502    -193.307    502    -193.256    502    -193.13    502 

    -203.948    500    -204.114    500    -204.047    500    -203.927    500 

    -192.587    504    -192.632    504    -192.616    504    -192.576    504 

    -191.827    498    -192.032    498    -191.962    498    -191.773    497 

    -194.707    500    -195.243    500    -195.057    500    -194.598    499 

    -195.202    498    -195.243    498    -195.227    498    -195.194    498 

    -195.921    496    -195.927    496    -195.921    496    -195.924    496 

    -195.392    498    -195.356    498    -195.365    498    -195.4    498 

    -192.545    498    -192.6    498    -192.58    498    -192.532    498 

    




and for 10 cells


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -139.783    11    -147.468    10    -143.755    11    -139.270    11 

    -137.910    10    -142.610    10    -140.283    10    -137.284    10 

    -142.794    10    -146.055    9    -144.733    9    -142.067    10 

    -133.883    9    -134.231    10    -134.058    10    -133.530    9 

    -140.499    10    -142.933    10    -141.646    10    -140.285    10 

    -130.732    10    -129.973    10    -130.279    10    -131.016    10 

    -133.887    10    -133.595    10    -133.623    10    -134.086    10 

    -135.397    9    -133.938    9    -134.267    9    -136.617    9 

    -139.761    10    -141.720    10    -140.478    10    -140.018    10 

    -142.768    10    -147.142    9    -145.116    9    -142.328    10 



Again the estimates are good for all models, and similar
log-likelihoods are obtained.  Reducing ϕ to 0.07, we obtain for
500 cells





  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -177.336    495    -177.035    495    -177.027    495    -177.141    495 

    -175.239    477    -174.969    479    -175.003    479    -175.107    478 

    -176.421    497    -176.359    498    -176.295    498    -176.324    497 

    -171.300    460    -171.109    466    -171.179    464    -171.251    462 

    -176.719    537    -178.073    538    -177.447    538    -176.984    538 

    -184.978    501    -187.879    489    -186.498    494    -185.504    498 

    -177.330    523    -178.051    524    -177.673    524    -177.434    523 

    -183.672    510    -185.782    501    -184.699    504    -183.973    507 

    -171.51    523    -171.445    532    -171.511    529    -171.532    526 

    -176.99    516    -175.459    518    -175.868    517    -176.388    516 
 



and for 10 cells


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -90.0015    12    -112.940    4    -96.5969    6    -74.8928    10 

    -90.1253    13    -112.852    4    -96.6292    7    -77.1458    11 

    -78.6269    8    -96.6019    4    -83.227    5    -63.8753    8 

    -88.6604    12    -108.518    4    -93.0871    7    -73.4382    10 

    -80.5898    11    -93.9211    3    -80.5046    4    -66.668    9 

    -86.3204    9    -107.53    5    -93.3741    7    -81.5557    9 

    -69.2068    8    -80.6773    2    -68.6397    3    -52.4727    7 

    -77.0180    9    -92.9245    3    -79.659    5    -61.3999    8 

    -81.5496    8    -108.38    4    -92.4763    6    -74.1947    9 

    -91.0582    11    -106.786    4    -92.8809    6    -78.1598    10 
 
  



We see that variability in the cell estimates is increasing, and that
for the 500-cell simulations the log-likelihoods are close amongst the
four models in each simulation. However for the 10-cell simulations we
see that the based model has higher maximized likelihoods in all
the simulations, indicating a better overall fit to the simulated peak
heights. We can also see that the lognormal and gamma models appear to
be underestimating the true cell-count value, whilst the normal and
models are less biased in their estimates, but the normal model
having a greater variability. The following boxplot of cell counts
estimated for 10 cells, and ϕ = 0.07, is based on 200
simulations, and confirms the small sample behaviour in the table
above:



  
    < g r a p h i c s >



 



  §.§.§ Simulations with stutter but no noise 


We now include stutter products in the simulations. For all
simulations we set the conditional probability of stutter on each
allele to be 0.004, and for forward and double stutter conditional
probabilities are all set to 0.001. We do not include background
noise, but we set the analytic threshold to be 30 .

For ϕ = 1 and 500 cells, we obtain the following estimates from
10 simulations. We see that all models give similar cell count
estimates and log likelihoods. However we see that there is much more
variation amongst the maximized likelihoods between the four models in
each simulation.


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -398.925    500    -397.548    500    -397.95    500    -400.347    501 

    -401.835    500    -400.62    500    -400.966    500    -400.65    500 

    -390.493    501    -389.57    501    -389.845    501    -396.683    502 

    -401.541    500    -401.002    500    -401.14    500    -403.778    500 

    -399.756    500    -399.765    500    -399.731    500    -401.508    500 

    -390.409    499    -390.135    499    -390.196    499    -397.262    500 

    -403.813    500    -402.911    500    -403.161    500    -400.463    501 

    -394.292    499    -393.968    499    -394.046    499    -399.208    500 

    -393.24    498    -392.758    498    -392.896    498    -396.412    499 

    -392.309    500    -392.337    500    -392.311    500    -401.011    500 

    




Reducing the number of cells to 200 we obtain the following table. We
see that all models are making accurate estimates of the numbers of
cells, but we also see that the model likelihoods are, with one
simulation exception, lower than the moment based models.


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -344.81    200    -343.656    200    -344.009    200    -354.253    200 

    -349.42    199    -346.64    199    -347.456    199    -353.245    200 

    -350.053    199    -348.586    199    -348.993    199    -354.592    199 

    -365.401    200    -364.874    200    -364.998    200    -372.019    200 

    -343.929    200    -342.793    200    -343.114    200    -351.199    200 

    -369.769    200    -366.598    200    -367.482    200    -366.688    200 

    -342.514    200    -341.61    200    -341.877    200    -350.469    200 

    -335.855    201    -335.204    201    -335.411    201    -347.73    201 

    -344.183    200    -343.186    200    -343.482    200    -349.88    200 

    -358.427    199    -357.657    199    -357.856    199    -363.16    200 

    



Reducing the number of cells down to 100 we obtain the following
table.  Cell estimates are again good for all models, but we see that
divergence between the log-likelihood values and the moment
models is becoming more pronounced.



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -233.411    100    -233.38    100    -233.384    100    -240.193    100 

    -236.636    101    -235.228    101    -235.645    101    -235.647    101 

    -237.502    100    -236.287    100    -236.657    100    -240.525    100 

    -234.364    99    -234.535    99    -234.442    99    -236.879    99 

    -223.076    100    -223.137    100    -223.106    100    -229.342    100 

    -228.312    101    -227.795    101    -227.945    101    -231.916    101 

    -230.257    99    -229.443    99    -229.69    99    -232.93    99 

    -228.276    100    -227.977    100    -228.071    100    -234.852    100 

    -224.315    100    -223.785    100    -223.957    100    -231.222    100 

    -221.219    100    -221.08    100    -221.128    100    -229.573    100 

    



Reducing the number of cells down to 10 we obtain the following table,
in which all models correctly estimate the number of cells in each
simulations, but the divergence between the log-likelihood values
and the moment models is still apparent.


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -108.134    10    -108.002    10    -108.038    10    -115.534    10 

    -111.576    10    -111.361    10    -111.413    10    -119.017    10 

    -110.092    10    -110.252    10    -110.185    10    -117.129    10 

    -111.932    10    -111.352    10    -111.528    10    -119.666    10 

    -110.431    10    -111.194    10    -110.915    10    -117.359    10 

    -107.527    10    -108.362    10    -108.054    10    -114.481    10 

    -106.965    10    -106.724    10    -106.799    10    -114.368    10 

    -113.55    10    -113.298    10    -113.352    10    -121.412    10 

    -109.083    10    -108.913    10    -108.956    10    -116.539    10 

    -116.432    10    -117.402    10    -117.013    10    -123.514    10 

    




We now set ϕ = 0.2.  For 500 cells, we obtain the following
table, in which we are now seeing the effects of the pre-sampling
variability in the wider range of cell estimates in all models, which
are giving similar estimates.



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -279.819    496    -280.204    497    -280.053    497    -284.162    498 

    -313.331    497    -307.348    493    -308.797    494    -302.683    495 

    -279.83    496    -279.6    498    -279.678    497    -283.113    498 

    -279.828    504    -279.664    506    -279.717    505    -283.564    506 

    -288.832    496    -288.937    496    -288.811    496    -290.632    497 

    -302.379    497    -298.563    495    -299.36    495    -297.065    496 

    -286.874    504    -287.079    504    -286.944    504    -289.078    505 

    -286.703    507    -286.583    507    -286.567    507    -289.549    509 

    -301.076    502    -301.912    500    -301.295    500    -299.415    500 

    -292.959    500    -293.018    499    -292.885    500    -294.815    501 

    




Reducing the number of cells to 200, and keeping ϕ = 0.2 we
obtain



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -227.224    199    -228.354    199    -227.737    199    -228.009    198 

    -237.356    196    -233.14    195    -234.183    195    -233.538    195 

    -229.561    205    -231.812    206    -230.937    206    -232.314    207 

    -233.089    196    -237.701    193    -235.661    194    -234.839    195 

    -223.573    197    -222.67    199    -222.915    198    -224.836    198 

    -233.505    205    -233.881    205    -233.493    205    -234.295    205 

    -232.002    198    -230.505    197    -230.724    197    -231.164    196 

    -225.549    195    -225.824    196    -225.692    195    -227.969    196 

    -216.563    192    -218.336    194    -217.713    193    -219.421    193 

    -220.548    199    -221.047    200    -220.84    199    -222.54    200 

    



Reducing the number of cells to 100, and keeping ϕ = 0.2 we
obtain similar behaviour,




  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -161.866    91    -163.224    92    -162.538    92    -162.197    91 

    -174.273    102    -174.564    99    -174.032    100    -173.428    100 

    -164.162    99    -163.554    100    -163.703    100    -164.374    99 

    -166.503    100    -168.923    100    -167.691    101    -167.244    100 

    -165.579    96    -165.367    96    -165.226    96    -165.526    96 

    -161.76    99    -161.274    101    -161.398    100    -161.841    100 

    -169.946    106    -171.175    105    -170.44    105    -170.32    106 

    -164.081    98    -164.207    99    -164.096    99    -164.48    98 

    -165.631    98    -167.698    97    -166.658    98    -166.109    98 

    -174.81    103    -173.074    101    -173.319    102    -174.481    102 

    



However as we go below 40 cells the behaviours of the and moment
based models start to change.  For 40 cells we obtain the following
table, in which all the models are in broad agreement:



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -147.227    38    -151.059    38    -149.056    39    -147.405    39 

    -149.565    44    -149.253    44    -149.071    44    -148.65    46 

    -144.248    39    -145.536    40    -145.022    40    -144.202    41 

    -149.696    40    -159.711    38    -154.874    39    -150.288    41 

    -148.827    41    -148.863    41    -148.599    41    -148.133    42 

    -149.489    41    -149.389    41    -149.102    41    -148.598    42 

    -151.896    41    -153.712    40    -152.28    40    -150.926    41 

    -142.965    33    -143.202    35    -142.985    35    -143.77    37 

    -145.254    36    -149.028    36    -146.878    36    -145.521    37 

    -153.346    37    -149.71    35    -150    36    -151.28    37 
 
    



However the further reduction to 35 cells produces the following table
in which we start to see the model overestimating the number of
cells in a biased manner, although the likelihoods of the four models
are in broad agreement.


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -144.129    36    -148.72    36    -146.337    37    -144.584    37 

    -146.785    36    -149.664    36    -147.998    37    -146.792    37 

    -150.221    33    -151.583    32    -149.991    33    -151.201    37 

    -142.658    34    -143.072    36    -142.915    36    -143.154    37 

    -152.339    33    -153.707    30    -152.106    31    -153.466    35 

    -147.424    37    -144.977    38    -145.683    38    -146.49    39 

    -152.856    36    -156.567    34    -154.021    35    -152.478    37 

    -152.798    33    -154.999    28    -152.755    30    -155.032    35 

    -146.274    34    -150.299    34    -148.062    34    -147.571    37 

    -144.768    36    -145.235    37    -144.882    37    -144.767    37 

    


Reducing the number of cells to 30 yields divergence in the likelihood
estimates too, as in the following table:


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -140.094    28    -142.829    29    -141.302    29    -145.536    34 

    -140.487    29    -143.953    29    -142.178    29    -146.174    34 

    -137.771    27    -138.535    29    -138.321    28    -145.971    33 

    -142.205    30    -145.612    28    -143.452    29    -146.512    34 

    -141.715    30    -146.407    29    -143.991    30    -146.539    34 

    -147.524    30    -146.787    29    -146.315    29    -151.487    34 

    -145.765    30    -146.033    30    -145.324    30    -149.674    34 

    -139.777    30    -143.263    30    -141.447    30    -144.618    34 

    -146.34    32    -153.725    28    -149.89    29    -150.328    35 

    -148.841    31    -147.083    30    -146.882    31    -150.594    35 

    



Reducing the number of cells to 10 we obtain a stark difference
between the and moment based models, as shown in the following
table:



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -73.3783    10    -76.5299    10    -75.1702    10    -138.845    23 

    -69.2668    10    -69.0933    10    -68.4735    10    -129.504    23 

    -66.276    10    -63.7138    10    -64.0525    10    -130.571    23 

    -83.0595    11    -84.5144    13    -84.3105    12    -138.676    24 

    -67.5142    9    -68.668    10    -68.2279    9    -137.273    3 

    -81.0544    12    -81.9969    12    -81.2933    12    -132.806    24 

    -80.4326    11    -79.8365    11    -79.5108    11    -137.78    23 

    -66.6831    9    -67.4642    11    -67.4318    10    -130.458    23 

    -76.19    11    -78.2126    11    -77.0102    11    -135.338    23 

    -65.9317    9    -65.3378    10    -65.1538    10    -130.649    23 

    



It appears that the moment based models are giving close to unbiased
estimates for the number of cells, whereas the model is giving
quite biased estimates. The problem with the model can be traced
to factorisation approximation of (<ref>).  Recall that
this is a double approximation here: one in which the joint is
replaced with a factorised product, given in (<ref>),
and the second factorisation approximation made in going to
(<ref>). Now consider an allele having k genomes. The
probability that it completely drops out will be
(1-ϕ)^k. However, if the allele completely drops out then it
cannot produce any stutter products, hence the joint probability of
dropout for that allele and its stutter products will be
(1-ϕ)^k. However the factorisation approximation
(<ref>) means that the factor (1-ϕ)^k is counted
for the allele and each of it stutter product, hence will give a total
joint probability of dropout of (1-ϕ)^4k. For large ϕ and
k, both (1-ϕ)^k and (1-ϕ)^4k will be very close to zero,
and so the factorisation approximation will be good, but if k and/or
ϕ are small, then the terms will have a large effect on the
likelihood calculations, and their differences will also grow.  It
appears from the results in the tables above that the transition point
where the approximation starts becoming bad is when the number of
cells is around 35-40: with ϕ=0.2 this corresponds to around 7-8
amplify-able cells, approximately 45-50pg by weight, of DNA
subject to .

It is curious that the moment models appear to give approximately
unbiased estimates for the number of cells.  If we compare to the
10-cell, ϕ=1 table on page:cell10tab, in which the
lognormal and gamma models gave underestimates for the number of
cells, it appears that the factorisation approximation has had a
correcting effect on these models. The normal model appears to have
less variability in its log-likelihood values when compared to to the
lognormal and gamma models.

Now ϕ=0.2 is a high value for forensic samples. Recall that
ϕ is the product of the extraction efficiency and the fraction of
extracted DNA that is taken for amplification.  A more typical
extraction efficiency value would be around 0.2, and similarly the
fraction of aliquot taken for amplification would also be around 0.2,
so that ϕ values of around ϕ = 0.04 would be more realistic
for forensic casework <cit.>. Hence, if the transition
for the model happens when expected number of amplify-able cells
is around 7-8, this corresponds to around 200 cells in the sample.

The following table shows estimated cell counts and log-likelihoods
for 200 cells and ϕ = 0.04, corresponding to around 8 cells of
DNA. We see that all the models are quite comparable in their
estimates of cell counts and their likelihoods values.



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -153.806    218    -158.376    212    -155.638    217    -154.016    222 

    -158.168    213    -160.208    199    -158.317    206    -157.783    208 

    -156.432    212    -158.452    204    -156.755    208    -156.127    208 

    -154.188    195    -155.784    189    -154.183    193    -153.931    195 

    -153.257    201    -152.824    199    -152.443    201    -152.747    202 

    -145.644    177    -143.678    190    -144.625    187    -146.097    190 

    -157.511    206    -157.945    190    -156.17    196    -155.666    202 

    -149.454    197    -149.449    203    -149.159    202    -148.911    202 

    -146.514    183    -152.516    182    -149.646    185    -148.237    192 

    -156.928    227    -157.76    220    -156.721    224    -155.927    224 

    




The following boxplots for estimated cell counts and log-likelihoods,
based on 200 simulations confirms this:



  .5
    
    < g r a p h i c s >

  
  .5
    
    < g r a p h i c s >

  



However, reducing the number of cells to 150 we obtain the following
table


  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -144.025    136    -142.352    143    -142.807    142    -151.579    170 

    -137.712    141    -144.199    131    -140.921    136    -146.336    170 

    -140.71    141    -142.032    147    -141.36    147    -147.833    173 

    -150.072    171    -151.597    164    -150.225    167    -151.937    188 

    -145.085    138    -142.67    143    -143.295    143    -151.939    171 

    -140.398    135    -142.12    134    -140.728    136    -148.936    168 

    -144.896    146    -146.134    145    -144.956    147    -151.04    175 

    -142.788    140    -145.464    127    -143.368    133    -151.46    169 

    -137.54    131    -141.125    134    -139.346    135    -147.239    167 

    -148.65    156    -146.808    158    -147.003    158    -151.945    181 

    


and these corresponding plots based on 200 simulations:


  .5
    
    < g r a p h i c s >

  
  .5
    
    < g r a p h i c s >

  



Reducing the number of cells to just 30, (so giving an average of just
1.2 cells worth of DNA for ) we obtain the following:



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -26.6781    26    -26.6658    27    -26.226    24    -51.4625    4 

    -31.7871    30    -31.8301    32    -31.5047    29    -64.4055    6 

    -29.1594    27    -29.4184    33    -29.4628    29    -61.8782    5 

    -35.5867    28    -37.4363    39    -37.6121    34    -72.5157    6 

    -25.1085    25    -24.6105    28    -24.5907    24    -51.6145    4 

    -20.4599    24    -20.3479    25    -20.191    22    -44.142    4 

    -27.5168    25    -28.6467    30    -28.5487    26    -53.8391    4 

    -35.9236    31    -36.4142    35    -36.2378    32    -74.6161    7 

    -17.6097    21    -17.9368    25    -18.0696    21    -38.5006    3 

    -38.5072    33    -37.1138    34    -36.9363    31    -75.1298    7 

    



and these corresponding plots:


  .5
    
    < g r a p h i c s >

  
  .5
    
    < g r a p h i c s >

  



However, if we look at the simulated data for these latter
simulations, we find that they have complete drop-out on most of the
markers, and just one peak on each of the other markers. For such
samples it is not necessary to include double and forward stutters in
the model.  If we remove these model components, so that we include
only the target and single back-stutter allele distributions, we
obtain the following fits for the models (using the same simulated
data).




  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -15.3989    21    -14.676    18    -14.4406    16    -20.3993    11 

    -26.7267    26    -26.7212    27    -26.2684    24    -32.4348    25 

    -31.7885    30    -31.8266    32    -31.4997    29    -37.5075    31 

    -29.1259    28    -29.3837    33    -29.4362    29    -35.7485    31 

    -35.5782    28    -37.4145    39    -37.5979    34    -42.0806    33 

    -25.0822    25    -24.588    28    -24.5746    24    -31.628    26 

    -20.4623    24    -20.3585    25    -20.1993    22    -27.3342    23 

    -27.4902    25    -28.6197    30    -28.5286    26    -33.9025    25 

    -35.9469    31    -36.4683    35    -36.2816    32    -41.3409    35 

    -17.6271    21    -17.9428    25    -18.0737    21    -24.178    19 

    




We see that the estimate are now much better. The lower
log-likelihoods are explained by the inclusion of modelling
stutters. However the few observed peaks heights in the simulations
are within the range 30-100, so that any stutter peaks would not be
observed with the threshold 30. So, if we also remove the modelling of
stutter, then we obtain the following fits, in which we see that the
log-likelihoods of all four models are now comparable.




  
    2|c|Normal     2c|Logormal     2c|Gamma     2c| 

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells 
 -15.366    21    -14.6613    18    -14.4301    16    -14.3498    15 

    -26.6653    27    -26.6806    27    -26.2362    24    -26.1813    22 

    -31.7059    30    -31.7632    32    -31.4469    29    -31.785    27 

    -29.2184    28    -29.4225    33    -29.4644    29    -30.134    27 

    -35.5742    29    -37.4283    39    -37.6174    34    -36.6476    29 

    -25.388    26    -24.8448    28    -24.7658    24    -25.9578    23 

    -20.4145    24    -20.3252    25    -20.1738    22    -21.0171    21 

    -27.437    25    -28.5662    30    -28.4889    26    -27.6955    23 

    -35.856    32    -36.3818    36    -36.2137    32    -36.0407    31 

    -17.5942    21    -17.9095    25    -18.0512    21    -17.7821    19 

    








  §.§.§ Simulations with degradation, but no baseline noise


In the simulations above, we have simulated samples with a degradation
parameter equal to zero, and have used that value in estimating the
cell counts. We now consider simulations in which the degradation
parameter is non zero, and is estimated from the data along with the
cell counts.

Our first simulation set will have ϕ=0.04, and degradation
parameter δ = 0.01. This means that for an allele of size l
in base-pairs, that its binomial selection probability if
ϕexp(-δ l).  So an an allele of size 100 will have a
selection probability of around 0.0147, whilst a longer allele of
length 400 will have a selection probability of around 0.0007.

We take the initial number of cells to be 1000. Including all possible
forms of stutter in fitting the model we obtain the following
estimates from 10 simulations (the degradation for only the FFT model
is shown);



  
    2|c|Normal     2c|Logormal     2c|Gamma     2c      δ̂_FFT

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells     
 -119.379    1055    -120.355    1003    -119.455    1029    -147.95    1172    0.00995539 

    -110.236    1002    -111.833    1038    -111.143    1030    -140.249    1146    0.0100925 

    -118.145    1218    -116.529    1127    -115.724    1154    -147.429    1320    0.0111535 

    -123.923    1033    -123.257    1113    -123.683    1092    -149.485    1172    0.00954747 

    -113.382    1020    -111.865    1012    -111.872    1017    -142.634    1147    0.0101153 

    -128.499    1108    -127.712    1028    -126.872    1059    -153.195    1204    0.00958187 

    -133.743    1002    -133.587    954    -132.82    965    -161.607    1098    0.00973301 

    -116.822    1038    -117.803    1036    -116.58    1037    -146.971    1165    0.0103188 

    -124.35    1025    -122.464    1058    -122.859    1049    -149.219    1145    0.0093805 

    -117.821    1082    -114.408    1066    -114.787    1066    -148.147    1207    0.0108218 

    



Give the quite high level of degradation we would not expect to see
double or forward stutters, even amongst the smaller alleles.  So
removing these from the model, but retaining the stutter model, we
obtain (for the same simulated data) much better concordance between
the model and the moment-based models.




  
    2|c|Normal     2c|Logormal     2c|Gamma     2c      1c|δ̂_FFT

    LL_max     cells    LL_max     cells    LL_max     cells    LL_max     cells     
 -119.262    1057    -120.176    1007    -119.301    1033    -120.526    1058    0.00995539 

    -110.268    1004    -111.917    1040    -111.211    1032    -111.825    1027    0.0100925 

    -118.006    1220    -116.363    1132    -115.572    1159    -117.445    1190    0.0111535 

    -123.927    1035    -123.266    1116    -123.691    1095    -124.538    1070    0.00954747 

    -113.252    1023    -111.715    1016    -111.741    1021    -113.954    1027    0.0101153 

    -128.584    1111    -127.824    1030    -126.974    1061    -128.595    1094    0.00958187 

    -133.737    1005    -133.594    957    -132.831    968    -133.919    984    0.00973301 

    -116.784    1040    -117.731    1039    -116.525    1040    -117.751    1046    0.0103188 

    -124.403    1028    -122.545    1060    -122.924    1052    -124.78    1043    0.0093805 

    -117.763    1084    -114.387    1069    -114.757    1070    -117.423    1076    0.0108218 

    



The following plot shows the peak height data used in the last
simulation of the previous table. We see that there is complete
allelic drop-out on the loci CSF1PO and D18S51.



  
    < g r a p h i c s >



 
The following plot shows the peak height probability distribution
obtained from the model, for allele 25 of the locus FGA, in which
the distribution is conditional on the values all of the other peak
heights and the profile of the contributor. The smaller vertical line
at (30,0) locates the analytic threshold, the larger vertical line
locates the observed peak. The complete dropout probability value
point at (0, 0.0569), is off the scale of the plot. The plot
suggests that, most likely, either 3 or 4 intact genomic strands were
randomly selected for the simulated amplification.



  
    < g r a p h i c s >



 



  §.§.§ Summary for single contributor simulations

For moderate to large numbers of cells, all three moment models and
the based model yield comparable estimates for estimated cell
counts and maximized log-likelihood values. However when the expected
number of cells drops below 8 or so, the model exhibits bias in
overestimating the number of cells, and produces more extreme (lower)
likelihoods.  Overall the lognormal and gamma models appear to give
the better estimates. However, for such low amounts it seems sensible
to simplify the model by omitting the forward and double backward
stutter model components; doing so it appears that the model is
concordant with the moment based models.



 §.§ Two person simulations



We now look at a few two-person simulations, concentrating on
low-template scenarios with various relative initial amounts of DNA.



  §.§.§  Two person mixtures


In the first simulation we have no degradation, that is δ = 0,
set ϕ = 0.04, and we simulate from contributors
 (C_1=200 cells) and 
(C_2=200 cells).  With ϕ=0.04 this corresponds to an average of
16 amplify-able cells worth of DNA, around 150pg.  Using all four
stutter distributions, we obtain this set of simulated fits, assuming
that both contributor profiles are known. All models are giving
similar maximized likelihood and cell estimates.



  
    3|c|Normal     3c|Logormal     3c|Gamma     3c|  

    LL_max     Ĉ_1     Ĉ_2     LL_max      Ĉ_1     Ĉ_2    LL_max      Ĉ_1     Ĉ_2    LL_max      Ĉ_1     Ĉ_2  
 -242.429    193    194    -241.955    200    201    -242.155    198    200    -242.31    200    200 

    -244.799    202    196    -241.404    210    198    -242.383    208    198    -243.576    208    200 

    -248.631    184    181    -248.005    182    180    -247.48    184    181    -249.832    201    196 

    -246.318    175    207    -244.447    187    206    -244.682    182    206    -246.326    194    206 

    -249.793    207    224    -247.785    207    226    -247.889    208    226    -248.492    208    224 

    -240.07    180    214    -239.761    185    225    -239.775    185    222    -240.727    190    222 

    -253.473    204    212    -251.879    200    209    -251.465    203    210    -251.749    204    213 

    -252.464    197    207    -251.199    192    206    -250.913    194    208    -251.412    195    208 

    -246.45    198    214    -252.804    210    197    -249.38    206    205    -247.265    205    206 

    -245.234    191    209    -243.718    200    209    -243.995    198    210    -244.59    200    212 

    





Repeating the simulation but with C_1= 50 and C_2=200 cells, (a
1:4 mixture), we see that the three moment models are in broad
agreement, but that the FFT model is diverging in both the cell
estimates for  and the maximized likelihoods,
which are much lower than the moment models in comparison.



  
    3|c|Normal     3c|Logormal     3c|Gamma     3c|  

    LL_max     Ĉ_1     Ĉ_2     LL_max      Ĉ_1     Ĉ_2    LL_max      Ĉ_1     Ĉ_2    LL_max      Ĉ_1     Ĉ_2  
 -181.783    41    239    -180.183    40    237    -179.908    39    238    -215.229    9    251 

    -193.674    50    193    -192.84    47    193    -191.837    46    195    -239.032    16    204 

    -186.171    42    209    -190.008    51    212    -188.724    47    213    -232.002    13    224 

    -193.501    50    199    -190.686    54    200    -191.185    52    200    -234.826    111    189 

    -188.659    48    200    -191.771    46    204    -189.964    47    204    -235.36    15    204 

    -201.354    56    202    -200.354    60    196    -199.818    57    200    -240.015    115    192 

    -172.509    34    201    -173.196    41    215    -173.62    38    212    -205.595    9    208 

    -177.388    41    189    -176.959    42    196    -176.441    40    195    -213.534    11    204 

    -190.686    44    201    -194.281    48    188    -191.67    45    195    -233.393    13    204 

    -194.286    52    196    -200.514    54    192    -197.259    53    195    -237.177    112    187 

    



However if we take out the double and forward stutter models, we see
that balance is restored:



  
    3|c|Normal     3c|Logormal     3c|Gamma     3c|  

    LL_max     Ĉ_1     Ĉ_2     LL_max      Ĉ_1     Ĉ_2    LL_max      Ĉ_1     Ĉ_2    LL_max      Ĉ_1     Ĉ_2  
 -181.742    42    240    -180.308    42    238    -179.996    41    239    -183.262    45    237 

    -193.528    51    194    -192.87    49    193    -191.85    47    195    -195.064    52    193 

    -186.354    43    209    -190.206    52    213    -188.941    49    213    -191.354    50    212 

    -193.559    51    199    -190.807    55    201    -191.313    53    201    -194.136    55    200 

    -188.622    49    201    -191.907    48    204    -190.039    48    204    -191.171    53    201 

    -201.419    57    203    -200.684    61    196    -200.109    59    200    -202.325    60    202 

    -172.836    35    202    -173.856    43    215    -174.24    40    212    -176.951    39    207 

    -177.56    42    189    -177.187    44    196    -176.674    41    196    -179.332    44    192 

    -190.541    44    202    -194.199    49    189    -191.58    46    196    -192.827    49    199 

    -194.748    53    196    -200.858    56    192    -197.676    55    195    -197.303    57    196 

    




We now look at the data from the final simulation. We refit the set of
models under three more scenarios, making the following four
scenarios, in which the K1K2 scenario is that use in the previous
table.



  K1K2 scenario  genotype is treated as
  known, and  genotype is treated as known

  U1K2 scenario  genotype is treated as
  unknown, and  genotype is treated as known

  K1U2 scenario  genotype is treated as
  known, and  genotype is treated as unknown

  U1U2 scenario  genotype is treated as
  unknown, and  genotype is treated as unknown


The genotypes of the individuals, and the observed peaks using the
analytic threshold of 30 RFUs, are shown in tab:sim2data, and
the simulated is plotted in fig:sim2pepg.  From
tab:sim2data we can see that a number of alleles have dropped
out, and that for some loci there is complete drop-out of the
low-donor contributor.









The results of fitting each of the models to each of the four
scenarios is shown in tab:sim2results. We see that for each,
the models give similar results, even with the considerable allelic
drop-out that is taking place.  Note also that the table may be used
to find for prosecution to defence likelihood ratios.  Suppose that
the prosecution case is that  is a contributor,
and the defence hypothesis that  is not. Then the
log-likelihood ratio in favour of the prosecution hypothesis will be
( -194.748 - (-209.557))/log(10) = 6.43 Bans.
tab:sim2ratioresults shows various combinations of
log-likelihood ratios.  Notice that the values are higher for the last
two rows, which are hypotheses concerning the presence of the major
contributor, than equivalents hypothesis comparisons in the first two
rows concerning the presence of the minor contributor.  This is to
anticipated. What is perhaps surprising are the quite high values in
the first two rows, concerning the presence of the minor contributor,
given the very low template DNA from  the minor contributor.







§ APPLICATION TO SAMPLE DATA



This section describes the performance of the model above to the
publicly available DNA dataset from the PROVEDIt Initiative
<cit.>. This is a very large
dataset of laboratory controlled single source and mixed DNA samples,
with amplifications carried out for the Identifiler
Plus™, PowerPlex16H™ and
Globalfiler™ kits. From it we shall use the sets of
samples for the Identifiler Plus™ kit that were
amplified for 28 cycles. We shall use the processed data in the EXCEL
files available from the PROVEDIt Initiative website, rather than the
raw fsa files that are also available. The EXCEL files contain allele
calls (including O/L designations) in Genemapper output format, in
which an analytic threshold has been set to 1 . Although these
files contain artefacts such as split-peaks and dye-blobs, and the
PROVEDIt Initiative provides an EXCEL spreadsheet to help remove
these, they have been left in so that the data can be processed in
batch automatically. Hence our analyses presented here are using quite
noisy data.






 §.§ Calibration of model parameters to the data




The samples from the PROVEDIt Initiative were prepared from the
dilution of high-density extracted DNA. Hence a Poisson model for the
amount of DNA would be appropriate. However, we use a binomial model
as presented so far, as this is more appropriate for real life
samples.[Essentially we are approximating Poisson
  distributions by binomial distributions, usually the approximation
  is the other way around.] We assume an extraction efficiency of
ψ = 0.3, a sample volume of 25μL of which 10μL is used in
the amplification, hence π_f = 0.06.  We also used a drop-in rate
per locus of 0.021.

To estimate the noise distribution, the peak heights of all
off-ladder (O/L) designated alleles in the high DNA single
source samples were extracted and treated as empirical distributions
for each lane - each lane distribution is based on over 40,000
values. A cut-off of 100 s was used for the noise distribution.
fig:idnoise shows the noise distributions for each of the
four dyes, omitting the long tail up to 100 s.




To form the data files, we retained for each sample all of the
alleles that were called, that is, those that were given an allele
designation that was not given as O/L.


Amplification probabilities were assigned as follows. It is assumed
that all strand probabilities in any given locus are the same, that is
p_g = p_g_d = p_h = p_h_d = p_a = p_a_d = p.  From the
high-template DNA single source, non-degraded, samples, in each sample
the peak heights of called alleles were added together. These total
peak height values were scaled by the estimated amount of initial DNA
in each sample (this information is provided in the name of each
sample). The mean peak heights of every locus in the kit was then
found. It was found that TPOX had the largest mean value.  This locus
was - somewhat arbitrarily - assigned an amplification probability of
0.85 for all types of strand. Using the fact that the number of cycles
was 28, and that the mean number of amplicons from a single strand is
approximately (1+p)^28, amplification probabilities for the other
loci were found so that their mean values scaled correctly according
to their empirical mean when compared to the mean peak height of TPOX.

Conditional stutter probabilities were either assigned values of 0.001
for forward and double stutter and 0.004 for single stutter, except
for the loci in the allelic ladder of the kit where mean stutter
ratios were available from the manufacturer's literature,
<cit.> in which case somewhat crude estimates
were obtained from the plots of stutter ratio versus allele lengths to
obtain stutter probabilities that matched in mean; forward and double
stutter probabilities for the alleles were obtained from these by
division by 4.


A scale factor of two million was used for all lanes to convert the
number of tagged amplicons to an equivalent value.

We shall focus on the samples that were amplified for 5 seconds, and
for these use an analytic threshold of 15 s.


With all of these parameters set, and using the laboratory estimated
amount in each sample, maximum likelihood estimation was used to
estimate the degradation and initial cell counts of each hypothesised
contributor.

 



 §.§ Analysis of a  two person sample

From the two-person data the author extracted the profiles of 189 samples.
These two person mixtures were treated as 3 person mixtures in
analysis.  The profiles are discussed here, and in later sections use
the naming convention of the PROVEDIt Initiative datasets.


Sample
  

This two sample mixture was prepared in 1 ratio of 1:9, and subject to
10 seconds of sonification to simulate degradation. The target amount
of DNA was 0.15ng. With our use of ϕ =0.06, this corresponds to
approximately 150/0.06 = 2500pg of DNA initially, that is around 378
cells in total, with around 22 cells worth in the amplified
sample. This is a highly degraded very low template sample, we
therefore omit the double and forward stutter components when fitting
the models: the plot is shown in fig:buexample1;



We first analyse this sample under the hypothesis that both of the
actual contributors are contributors, and an additional untyped
contributor is assumed to be present as well. Maximized likelihoods
and estimated cell amounts are shown in tab:buexample1. The
degradation is estimated on a discretized scale, and in this case all
four models have estimated the same degradation value.


Of particular note is that all models have estimated the cell amount
for the extra untyped contributor to be 0 - which is in agreement with
the manner of the preparation of the two-person mixture. The forensic
science literature contain much discussion concerning how to estimate
the number of contributors to a mixture(see for example
<cit.>).  In the models presented here the maximum
likelihood estimates can return 0 for cell counts, because the cell
amounts beign estimated are discrete integers. This is in marked contrast to other
continuous peak height models that model contributor amounts or
relative amounts by continuous variables, for which the nesting of
models by adding extra hypothesised contributors will lead to
increasing maximized likelihoods. <cit.>




<cit.> introduced
various statistically based diagnostics to check how well a model fits
the profile data. One such diagnostic is a QQ-plot in which for each
allelic peak h_a in the set H of peaks observed above the
threshold T the quantities
P( H< h_a H > T, h_b∈ H, b a) are calculated. If the
models is `true', then these values should follow a uniform
distribution, so that a plot of these sorted values should follow a
straight line when plotted against the quantiles of the uniform
distribution. fig:buexample1qqplot shows the QQ-plot for the
model, which is the model with the highest likelihood in
tab:buexample1: we see that the fit is very good.




As an example of the predictive distribution of a peak, we take the
allele 13 of locus D13S317. For that locus, three peaks were observed
above the threshold of 15 s: allele 11 with a peak height of
117, allele 12 with a height of 177 and allele 13 with a height of 27.
Contributor  has genotype (12/13) on locus
D13S317, and  has genotype (11/12). Thus the 13
peak can arise from the minor contributor 
(recall that we are not including forward stutter in this analysis) or
drop-in.  From the estimates in tab:buexample1, there are an estimated 71 initial
genomic strands of this allele type, sampled intact with probability
0.06 for amplification, hence an estimated mean of around 4.26
amplify-able strands .  fig:buexample1dist shows the
predictive distribution for this allele conditional on the observed
peak heights and estimated cell amounts. The first peak on the left is
the baseline-noise distribution for VIC lane (compare it to
fig:idnoise).  The plot is giving a clear signal that just
one amply-able strand of this allele was in the aliquot minitube when
the was carried out. For comparison, the distribution obtained
from using the Gamma model is also shown in
fig:buexample1dist.





tab:buexample1w3u shows estimates for each model analyzed
under the assumption that the mixture is of three untyped persons. The
gamma model has the highest likelihood, and now the normal model is
giving a non-zero cell count estimate for all three persons.
Distribution plots for allele 13 of locus D13S317 are given in
fig:buexample1dist3u.








 §.§ Analysis of a three person sample


The sample

has contributor amounts in the proportions of 1:9:9, with a medium
level of degradation cause by rDNase I enzyme added to the
extract. Analyzing as a four person mixture, assuming as knowns the
genotypes of the three contributors and one extra untyped individual,
and ignoring forward and double stutter, we obtain the estimates in
tab:buexample3w3k.  The
model correctly estimate zero cells for the untyped contributor,
and the model gives the highest likelihood.  The estimated
mixture proportions are in line with the experimentally prepared
ratios for all four models.  The QQ-plot for the model is shown
in fig:buexample2fftqqplot.




If we include forward and double stutter in the models, we obtain the
estimates obtained in tab:buexample3w3kfd. We see that the
model has become out-of-line with the other models in having a
much lower likelihood, and much lower estimates for the number of
cells from the minor contributor.  The normal model is now the one
having the highest likelihood: QQ-plots for the Normal and model
are shown in fig:buexample2normalqqplot.







§ A MODIFICATION OF THE MODEL TO CORRECT FOR CORRELATIONS


We have seen that in some cases we can improve the model by
removing the forward and double back-stutters; however this will not
be appropriate for all cases.  Here we introduce an adjustment to the
model that puts back some of the correlation that was discarded
in the factorization approximations of (<ref>) and
(<ref>); we denote the new model by .

Recall that the point of making those approximations was purely to
make the computations tractable, by enabling the peak height
likelihood to factorize into a product of convolutions of univariate
probability distributions of individual alleles. However, as pointed
out earlier, a side effect of this approximation is that the complete
dropout probability for an allele is included as a factor multiple
times in the likelihood. The modification here aims to reduce such
multiple factors.

For any given allele a, there will up to 6 distributions to convolve
in finding its peak height likelihood contribution:



  * The base-line noise distribution

  * The drop-in distribution

  * The peak height distribution arising from genomes of type a

  * The peak height distribution arising from stutter from genomes
  of type one repeat larger than a

  * The peak height distribution arising from double stutter from
  genomes of type two repeats larger than a

  * The peak height distribution arising from forward stutter from
  genomes of type one repeat smaller than a


In the new approximation we always retain the first three
distributions. Now consider the contribution from say stutter from one
repeat higher. Suppose that no peak above threshold is observed at the
one-repeat higher position. This means that any genomic material of
allelic type one-repeat higher than a has dropped out. In this case,
we can therefore omit the marginal stutter peak height distribution in
the convolution as it is extremely likely that any stutter product
would be too low to be observed.  However, given that stutter peaks
are typically around 5-15% of the target peak producing the stutter,
this means that even if a peak is observed at one-repeat higher than
a, but is below a small multiple of the analytic threshold then we
can assume that any stutter contribution it makes to the allele a
will not be seen.  A similar argument applies also to the double
stutter and forward stutter distributions—more so as these tend to
have smaller stutter ratios than single stutter peaks.

Hence the modification is as follows. In evaluating the peak height
likelihood for an allele a, when forming the marginal peak height
distribution:


  * Always include the base-line noise distribution.

  * Always include the drop-in distribution.

  * Always include the peak height distribution arising from genomes
  of type a.

  * Only include the peak height distribution arising from stutter
  from genomes of type one repeat larger than a if there is a peak
  above three times the threshold at that position.

  * Only include the peak height distribution arising from double
  stutter from genomes of type two repeats larger than a if there is
  a peak above three times the threshold at that position.

  * Only include the peak height distribution arising from forward
  stutter from genomes of type one repeat smaller than a if there is
  a peak above three times the threshold at that position.

  * Convolve all of the included distributions to obtain the
  marginal peak height distribution to use for finding the likelihood
  for the peak of the allele a.


A factor different to three for the multiple of the analytic threshold
could be used—further research is required to find the best factor
to use.



 §.§ Revisiting the two person simulation of
  sec:sim2pers


The following table shows the model fits obtained using the model for the simulation with C_1 = 50 cells and C_2 = 200, which
should be compared to the second table in sec:sim2pers: we
see that the model is now much more in line with the other
models both in its maximum likelihood values and the cell estimates of
the two contributors.





  
    3|c|Normal     3c|Logormal     3c|Gamma     3c|  

    LL_max     Ĉ_1     Ĉ_2     LL_max      Ĉ_1     Ĉ_2    LL_max      Ĉ_1     Ĉ_2    LL_max      Ĉ_1     Ĉ_2  
 -181.783    41    239    -180.183    40    237    -179.908    39    238    -180.334    38    237 

    -193.674    50    193    -192.84    47    193    -191.837    46    195    -191.74    47    199 

    -186.171    42    209    -190.008    51    212    -188.724    47    213    -188.69    44    215 

    -193.501    50    199    -190.686    54    200    -191.185    52    200    -193.364    52    202 

    -188.659    48    200    -191.771    46    204    -189.964    47    204    -188.564    47    204 

    -201.354    56    202    -200.354    60    196    -199.818    57    200    -205.226    61    202 

    -172.509    34    201    -173.196    41    215    -173.62    38    212    -172.733    34    210 

    -177.388    41    189    -176.959    42    196    -176.441    40    195    -176.003    39    195 

    -190.686    44    201    -194.281    48    188    -191.67    45    195    -190.566    44    200 

    -194.286    52    196    -200.514    54    192    -197.259    53    195    -195.177    51    199 

    





 §.§ Revisiting the three person mixture of
  sec:bu3person from the PROVEDit Initiative


Re-analyzing
,
with both forward and double stutter included using the model, we
obtain a maximized log-likelihood of -405.944809925, and estimated
cell counts of Ĉ_1 = 31, Ĉ_2 = 596, Ĉ_3 = 617
and Ĉ_4 = 0, less extreme than the estimates in
tab:buexample3w3kfd. The corresponding QQ-plot is little
changed from fig:buexample2normalqqplot and is omitted.




 §.§ Another 2-person mixture from the PROVEDit datatset


In all the examples seen so far, we have seen broad agreement between
the various models. We now look at an example from PROVEDit datatset
in which the models give quite different model fits. The dataset we
shall look at is . This
is a two-person mixture of 540pg total DNA template in the ratio of
1:9 from the two contributors. The U105 indicates that the sample was
subject to 105 minutes exposure to ultra-violet radiation (to degrade
the sample) before .  tab:bugood2epg shows the peak
heights of the sample, above the analytic threshold of 15 s, and
the profiles of the two contributors. A plot of the is shown in
fig:bugoodepg.

[htbp]
  

Peaks heights and the genotypes of constributors    .
 Locus             Allele/height      Allele/height      Allele/height      Allele/height     Allele/height      Allele/height      Allele/height      Allele/height 

 Amelogenin	   	X	/	X			   	X	/	Y	   	X	/	1686	   	Y	/	1049	   				   				   				   				   					
 
CSF1PO	   	10	/	12			   	10	/	12	   	9	/	33	   	10	/	848	   	10.2	/	17	   	11	/	53	   	12	/	867					   					

D13S317	   	12	/	13			   	11	/	12	   	10	/	67	   	11	/	1269	   	12	/	1125	   	
13	/	79					
 
D16S539	   	9	/	12			   	10	/	10	   	9	/	297	   	10	/	2682	   	11	/	22	   	12	/	183				
 
D18S51	   	17	/	18			   	18	/	19	   	15	/	15	   	16	/	18	   	17	/	130	   	18	/	834	   	19	/	633			

D19S433	   	14	/	14.2			   	13.2	/	14	   	12	/	18	   	12.2	/	92	   	13	/	88	   	13.2	/	1250	   	14	/	1295	   	14.2	/	141					

D21S11	   	30	/	31.2			   	27	/	29	   	26	/	29	   	27	/	616	   	28	/	61	   	29	/	786	   	30	/	122	   	31.2	/	100	   	35	/ 19	   36	/ 23	
 
D2S1338	   	22	/	23			   	19	/	22	   	18	/	92	   	19	/	1148	   	21	/	127	   	22	/	1230	   	23	/	109										
 
D3S1358	   	15	/	16			   	15	/	15	   	13	/	20	   	14	/	252	   	15	/	3613	   	16	/	266				

D5S818	   	12	/	13			   	11	/	12	   	10	/	49	   	11	/	848	   	12	/	989	   	13	/	124		

D7S820	   	11	/	11			   	8	/	9	   	8	/	523	   	9	/	426	   	11	/	138			

D8S1179	   	12	/	14			   	11	/	13	   	10	/	44	   	11	/	858	   	12	/	207	   	13	/	1309	   	14	/	204	
 
FGA	   	20	/	21			   	21	/	28	   	20	/	134	   	21	/	617	   	27	/	40	   	28	/	557	   	31.2	/	23				

TH01	   	9	/	9.3			   	8	/	9.3	   	7	/	90	   	8	/	1508	   	8.3	/	22	   	9	/	230	   	9.3	/	1896					

TPOX	   	8	/	10			   	8	/	11	   	7	/	23	   	8	/	1234	   	10	/	141	   	11	/	880						
 
vWA	   	16	/	17			   	13	/	17	   	12	/	37	   	13	/	1123	   	16	/	215	   	17	/	932					
 








We shall analyse this mixture under the assumption of two
contributors. If we assume the true genotypes of the contributors, and
include forward and double stutter, then obtain a maximized
log-likelihood of -756.241357805 using the normal model. Omitting the
forward and double stutter, the maximized log-likelihood increases
significantly to -651.9280003. Hence in all analyses we shall include
only stutter, and omit both forward and double stutter model
components.

tab:bu2goodresults shows the results of fitting  
models to four possible two-contributor scenarios.  The first line of
the table is the scenario in which we assume as known the profiles of
the (true) contributors  and
.  In the second line we replace
 with an untyped contributor labelled
.  In the third line we replace  with
an untyped contributor labelled .  In the last line we
treat both contributors as having unknown genotypes.

We see that cell estimates are broadly in line for all models except
the normal model, which has higher major contributor estimates, and in
agreement with the 1:9 preparation ratio of the sample. There is quite
a lot of variation between the models in their log-likelihood
estimates for each scenario, with the model outperforming the
moment based models, by a considerable margin. The log-likelihood
estimates are also out of line with the other models, in that for the
second and fourth scenarios in which the  has
been replaced by an unknown person, the likelihoods are lower than the
scenario in which both contributor genotype are assumed known,
indicating evidence against  being a contributor.
 



tab:bugoodratioresults shows various log-likelihood ratios
(expressed in bans) obtained by comparing pairs of scenarios from
tab:bu2goodresults. We see that all models are giving similar
log-likelihood ratios in favour of the major contributor, compared to
the maximum of 24.71 Bans (the inverse log profile probability of
, and the model gives much larger values
in favour of the minor contributor compared to the other models, with
the normal model suggesting very strong evidence against
 being a contributor (the inverse log profile
probability of  is 20.22 Bans).



fig:bugoodqqplots shows the QQ-plots for each of the models
under the scenario in which both contributors are known.


If we add a third un-profiled contributor, then all the models give
similar results. Here are the values under the scenario of the two
known contributors and a third untyped contributor.



  
    Model    LL_max     Ĉ_1     Ĉ_2     Ĉ_3 
 Normal    -489.24      153     1204     16 

    Lognormal     -500.16     153     1081     26 

    Gamma    -496.16     153     1121    27 

        -491.83     152     1140     17 
 



If we add another untyped contributor, so that we analyse the sample
assuming the profiles of the two (true) contributors and two untyped
contributors, we obtain the following estimates in which the model gives an estimate of no cells to the second untyped contributor,
unlike the moment based models:

  
    Model    LL_max     Ĉ_1     Ĉ_2 
        Ĉ_3     Ĉ_4 
 Normal    -487.71      150     1202     14    8

    Lognormal     -500.04     151     1085     17     13 

    Gamma    -495.48     150     1126    15   15 

        -491.83     152     1140     17    0
 







§ SUMMARY AND OUTLOOK

In this paper, a framework for modelling single source and multiple
donor forensic DNA samples has been elaborated for STR loci, based on
an idealisation of the steps for forming an from an initial DNA
sample. It was shown that multivariate probability generating
functions provide a succinct and efficient mathematical representation
of the steps in the process, and that their evaluation can be carried
out efficiently using . Factorization approximations to make the
mathematical models tractable were shown to introduce some biases for
low template samples, but these biases can be largely removed by a
modification, the model, which reintroduces, in a principled
manner, some of the correlation lost by factorization approximations
without affecting computational efficiency.

The examples treated in this paper have assumed that contributors are unrelated individuals. 
The framework may be extended to include relatedness of individuals, which might require that linkage between some 
 loci would need to be taken into account.

 An issue not explored in this paper is the robustness of the
 performance of the models to mis-specifications of model
 parameters. This could be readily explored via simulations in which,
 for example, data is simulated using one set of amplification,
 stutter and selection probabilities, and the simulated data is then
 fitted assuming different probabilities.  Given the reaonably good
 performance of the model on the experimental data from the PROVEDit
 initiative, for which reasonable values were assumed for the
 parameters but without any detailed knowledge of what the true
 experimental values were, one could anticipate that the model will
 prove to be reasonably robust to such parameter mis-specifications.


The current paper has focussed on obtained by capillary
electrophoresis, however it may not limited to this. Recently,
<cit.> showed that by substituting read coverages for peak
heights, autosomal SNP mixtures analyzed using massively parallel
sequencing can be interpreted by the open source software Euroformix
<cit.>, originally developed for s. By
making the same substitution, that is, by replacing the scale
factors by `read factors', the framework developed in this paper is
also applicable to such SNP mixtures. Moreover, because SNPs are
biallelic and do not stutter, the factorization approximations are no
longer required.  



Finally, other branching processes outside of the forensic applications focussed on in this paper can be treated by the methods
introduced in this paper.






PART: 
*Appendices


§ DERIVATION OF MOMENTS FOR AMPLICON MODEL


Let N denote the number of target amplicons, and M the number of
stutter amplicons, arising from the amplification of a single target
amplicon. We have that

    N n    = ∂ F_n(t,s)/∂ t|_t=1,s=1
     N(N-1) n    = ∂^2 F_n(t,s)/∂ t^2|_t=1,s=1
     M n    = ∂ F_n(t,s)/∂ s|_t=1,s=1
     N(N-1) n    = ∂^2 F_n(t,s)/∂ s^2|_t=1,s=1
     NM n    = ∂^2 F_n(t,s)/∂ s∂ t|_t=1,s=1

from which the variance N and M may be found, and hence the
correlation (N,M).  We take each in turn.



  §.§.§ N n:

Using the chain rule for differentiating, we have

    ∂ F_n(t,s)/∂ t =
    
    (1-p)∂ F_n-1(t,s)/∂ t +
    2p(1-ξ)F_n-1(t,s)∂ F_n-1(t,s)/∂ t + pξ∂ F_n-1(t,s)/∂ t G_n-1(s)

We now substitute t=1 and s=1, use the property that
F_n(1,1) = 1 for all n, to obtain a linear difference equation
that can be solved:


    N n     = (1-p)( N n-1) + 2p(1-ξ)( N n-1) + pξ ( N n-1) 
    
                    = (1 + p(1-ξ)) ( N n-1)
        = (1 + p(1-ξ))^n ( N 0)
       = (1 + p(1-ξ))^n


This could have been anticipated: a main peak allele is amplified with
probability p(1-ξ), hence starting with one allele and doing n
amplifications, the mean number is as shown, (1 + p(1-ξ))^n.



  §.§.§ N(N-1) n and N n:

These could also be written down directly based on the previous
result, however we go through the formal steps.

    N(N-1) n    = ∂^2 F_n(t,s)/∂ t^2|_t=1,s=1
    ∂^2 F_n(t,s)/∂ t^2   =
     
                                                 (1-p)∂^2 F_n-1(t,s)/∂ t^2  + 
                                                 2p(1-ξ)F_n-1(t,s)∂^2 F_n-1(t,s)/∂ t^2  +
       
                                                                                                                         2p(1-ξ)(∂ F_n-1(t,s)/∂ t)^2
                                                                                                                         + pξ∂^2 F_n-1(t,s)/∂ t^2 G_n-1(s)
    
       N(N-1) n    = (1 + p(1-ξ)) ( N(N-1) n-1) + 2p(1-ξ)( N n-1 )^2
    
       N^2 n    = (1 + p(1-ξ)) ( N^2 n-1) + 2p(1-ξ)( N n-1 )^2
    
                        = (1 + p(1-ξ)) ( N^2 n-1) + 2(p(1-ξ))(1+p(1-ξ))^2n-2
    
                        = (1 + p(1-ξ))^n ( N^2 0) + 2(p(1-ξ))∑_j=1^n (1+p(1-ξ))^2j-2(1 + p(1-ξ))^n-j
    
                        = (1 + p(1-ξ))^n + 2(p(1-ξ))∑_j=0^n-1 (1+p(1-ξ))^2j(1 + p(1-ξ))^n-1-j
    
                        = (1 + p(1-ξ))^n + 2(p(1-ξ))(1 + p(1-ξ))^n-1∑_j=0^n-1 (1+p(1-ξ))^j
    
                        = (1 + p(1-ξ))^n + 2(p(1-ξ))(1 + p(1-ξ))^n-1(1+p(1-ξ))^n-1/(1+p(1-ξ))-1
    
                        = (1 + p(1-ξ))^n [ 1 + 2(1+p(1-ξ))^n-1/1+p(1-ξ)]
    
                        = (1 + p(1-ξ))^n-1[ 2(1+p(1-ξ))^n +p(1-ξ)-1 ]
 
 
Hence the variance is

    N n    = ( N^2 n) - ( N n)^2
       =(1 + p(1-ξ))^n [ 1 + 2(1+p(1-ξ))^n-1/1+p(1-ξ)] -(1 + p(1-ξ))^2n
    
                   =(1 + p(1-ξ))^n [1 + 2(1+p(1-ξ))^n-1/1+p(1-ξ) - (1 + p(1-ξ))^n]
    
                   = 1-p(1-ξ)/1+p(1-ξ)(1 + p(1-ξ))^n [ (1+p(1-ξ))^n -1]



We also have

    N(N-1) n    = ( N n) + ( N n)^2 - ( N n) 
       = 1-p(1-ξ)/1+p(1-ξ)(1 + p(1-ξ))^n [ (1+p(1-ξ))^n -1]
                       +(1 + p(1-ξ))^2n - (1 + p(1-ξ))^n
    
                        = 2(1 + p(1-ξ))^n-1((1 + p(1-ξ))^n-1)




  §.§.§ M n:

This follows a similar pattern to the above for N:


    M n    = ∂ F_n(t,s)/∂ s|_t=1,s=1
    
      ∂ F_n(t,s)/∂ s   =
                                             (1-p)∂ F_n-1(t,s)/∂ s  + 2p(1-ξ)F_n-1(t,s)∂ F_n-1(t,s)/∂ s
       
          + pξ∂ F_n-1(t,s)/∂ s G_n-1(s)
          + pξ  F_n-1(t,s) ∂ G_n-1(s)/∂ s


Hence


    M n = (1 +p(1-ξ))( M n-1) + pξ∂ G_n-1(s)/∂ s|_s=1


Now 
    G_n(s) = (1-p)G_n-1(s) + pG_n-1^2(s),


which is the recursion relation for the generating function of a
single allele amplifying with probability p in the branching
process, hence


    M n = (1 +p(1-ξ))( M n-1) + pξ (1+p)^n-1


Now we start with no stutter amplicons, hence M 0 = 0. Thus


    M n    = (1 +p(1-ξ))( M n-1) + pξ (1+p)^n-1
       =  (1 +p(1-ξ))^2( M n-2) + (1 +p(1-ξ))pξ (1+p)^n-2+ pξ (1+p)^n-1
       = ⋮
       =  (1 +p(1-ξ))^n( M 0) + pξ∑_j=0^n-1(1 +p(1-ξ))^j(1+p)^n-j-1

We now use M 0 = 0, 1+p -(1+p(1-ξ) = pξ, and the
relation

    a^n-1 + a^n-2b + a^n-3b^2 + ⋯ ab^n-2 + b^n-1
    = a^n-b^n/a-b
 to deduce that


    M n = pξ(1+p)^n - (1+p(1-ξ))^n/(1+p) - (1+p(1-ξ)
    = (1+p)^n - (1+p(1-ξ))^n



When ξ=0 this vanishes as it should. Note also that
[N+M n] = (1+p)^n which is also as expected– it is as if the
stutter and the main target cannot be distinguished apart. We also
have


    M n / N+M n = 1 - ( 1 -pξ/1+p)^n


which → 1 and n →∞ for ξ > 0. Hence ultimately the
stutter allele peak will dominate the original allele peak with enough
cycles, a result found by <cit.>.



  §.§.§ M(M-1) n and M  n


Denote F_n;s(t,s) = ∂ F_n(t,s)/∂ s, and F_n;ss
the second derivative.

Hence

    M(M-1) n    = F_n;ss|_t=1,s=1
    
      F_n;ss   = (1-p) F_n-1;ss + 2p(1-ξ)F_n-1F_n-1;ss + 2p(1-ξ)F_n-1;s^2
       + pξ  F_n-1;ss G_n-1 +2 pξ F_n-1;sG_n-1;s + pξ F_n-1G_n-1;ss
    
       M(M-1) n    = (1 + p(1-ξ))( M(M-1) n-1) +
       (2p(1-ξ)F_n-1;s^2+2 pξ F_n-1;sG_n-1;s + pξ F_n-1G_n-1;ss)|_s=t=1
    
                        = (1 + p(1-ξ))( M(M-1) n-1) +
       
                       2p(1-ξ)((1+p)^n-1 - (1+p(1-ξ))^n-1)^2+
        
                       2 pξ((1+p)^n-1 - (1+p(1-ξ))^n-1)(1+p)^n-1 +
       
                       2pξ(1+p)^n-2[(1+p)^n-1 -1]


But we have that
M n = (1 +p(1-ξ))( M n-1) + pξ (1+p)^n-1, hence
adding this to both sides we obtain

    M^2  n    =  (1 +p(1-ξ))( M n-1) 
       +2p(1-ξ)((1+p)^n-1 - (1+p(1-ξ))^n-1)^2+
        
                     2 pξ((1+p)^n-1 - (1+p(1-ξ))^n-1)(1+p)^n-1 +
       
                     pξ(1+p)^n-2[2(1+p)^n-1 +p-1]

Expanding this out we have

    M^2  n    =  (1 +p(1-ξ))( M^2 n-1) 
       +2p(1-ξ)((1+p)^2(n-1) -2(1+p)^n-1(1+p(1-ξ))^n-1 + (1+p(1-ξ))^2(n-1))
       +
                     2 pξ( (1+p)^2(n-1) - (1+p)^n-1 (1+p(1-ξ))^n-1)
       +
                     pξ/1+p( 2(1+p)^2(n-1) -(1-p)(1+p)^n-1)

Collecting together similar terms, this reduces to

    M^2  n    =  (1 +p(1-ξ))( M^2 n-1) 
        +(2p + 2pξ/1+p)(1+p)^2(n-1)
        -
                     (4p(1-ξ) +2pξ))(1+p)^n-1(1+p(1-ξ))^n-1
       +2p(1-ξ)(1+p(1-ξ))^2(n-1)
       -p(1-p)ξ/1+p(1+p)^n-1

 

We may now solve this linear difference equation, noting that the
terms on the right are each of the form a^n-1 or a^2(n-1),
making the particular solution terms of the form (a^n-b^n)/(a-b) or
(a^2n-b^n)/(a^2-b), where in each case b = 1+p(1-ξ). We also
use that ( M^2 0) =0. We obtain
 
 
    M^2  n    = 
                      (2p + 2pξ/1+p)(1+p)^2n - (1+p(1-ξ))^n/(1+p)^2-(1+p(1-ξ))
        -  (4p(1-ξ) +2pξ))(1+p(1-ξ))^n-1(1+p)^n-1/1+p -1
        +2p(1-ξ)(1+p(1-ξ))^n-1(1+p(1-ξ))^n -1/(1+p(1-ξ)) -1
       -p(1-p)ξ/1+p(1+p)^n - (1+p(1-ξ))^n/(1+p)-(1+p(1-ξ))

 
 Now
 
    (1+p)-1 = p
    
       (1+p(1-ξ))-1 = p(1-ξ)
    
       (1+p) - (1+p(1-ξ)    = pξ
    
       (1+p)^2 - (1+p(1-ξ)   = 1 + 2p + p^2-1 - p +pξ = p(1+p+ξ)

 
 Hence simplifying further we obtain
 
    M^2  n    =  2(1+p)^2n - (1+p(1-ξ))^n/1+p
       -(4 - 2ξ)(1+p(1-ξ))^n-1 ((1+p)^n-1)
       + 2(1+p(1-ξ))^n-1((1+p(1-ξ))^n -1)
       -(1-p)/1+p ((1+p)^n - (1+p(1-ξ))^n)

 
 
 Now we had previously that


    M n =  (1+p)^n - (1 + p(1-ξ))^n


Squaring this and subtracting from E M^2 n gives the desired
variance M n. However there does not appear to be a nice
simple reduction in the result, so this is omitted.
 
 
For large n we may drop the terms in (1+p(1-ξ))^n and lower, and
we have approximately


    M^2 ≈ 2n (1+p)^2n-1 -
    (4-2ξ)(1+p(1-ξ))^n-1(1+p)^n-1 + 2(1+p(1-ξ))^2(n-1)

 
 
 
 


  §.§.§ NM n:


    NM n    = ∂^2 F_n(t,s)/∂ s∂ t|_t=1,s=1
    ∂^2 F_n(t,s)/∂ s∂ t   = 
                                                        (1-p)∂^2 F_n-1(t,s)/∂ s∂ t
                                                        +2p(1-ξ)F_n-1(t,s)∂^2 F_n-1(t,s)/∂ s∂ t
                                                        +pξ∂^2 F_n-1(t,s)/∂ s∂ tG_n-1(s)
      
       +2p(1-ξ)∂ F_n-1(t,s)/∂ s∂ F_n-1(t,s)/∂ t
       + pξ∂ F_n-1(t,s)/∂ t∂ G_n-1(s)/∂ s
     NM n   = (1+p(1-ξ))( NM n-1 )+2p(1-ξ)( N n-1)( M n-1)
                  +pξ ( N n-1)(1+p)^n-1
    
                    = (1+p(1-ξ))( NM n-1 )
    
                   +2p(1-ξ)(1+p(1-ξ))^n-1[ (1+p)^n-1-(1+p(1-ξ))^n-1]
      
       + pξ (1+p(1-ξ))^n-1(1+p)^n-1
    
      =   (1+p(1-ξ))( NM n-1 ) + p(2-ξ)(1+p(1-ξ))^n-1(1+p)^n-1
       
          -2p(1-ξ)(1+p(1-ξ))^n-1(1+p(1-ξ))^n-1



Using NM 0 = 0, the solution may be written as



    NM n    =  (2-ξ)p(1 + p(1-ξ))^n-1((1+p)^n-1)/(1+p-1)
       
                   -2p(1-ξ)(1 + p(1-ξ))^n-1(1 + p(1-ξ))^n-1/(1 + p(1-ξ) -1
       = (2-ξ)(1 + p(1-ξ))^n-1((1+p)^n-1) -2(1 + p(1-ξ))^n-1((1 + p(1-ξ))^n-1)
       = (1 + p(1-ξ))^n-1( (2-ξ)(1+p)^n-1) -2((1 + p(1-ξ))^n-1)
       = (1 + p(1-ξ))^n-1(ξ(1-(1+p)^n) + 2( (1+p)^n - (1 + p(1-ξ))^n)




  §.§.§ (M,N) 


The covariance (M,N) is given by

    (M,N)    =  NM n  - ( N n )( M n )
       = (1 + p(1-ξ))^n-1(ξ(1-(1+p)^n) + 2( (1+p)^n - (1 + p(1-ξ))^n)
       
                  - ((1+p)^n - (1 + p(1-ξ))^n))(1 + p(1-ξ))^n
    
                   = (1 + p(1-ξ))^n-1( (1-p(1-ξ))( (1+p)^n - (1 + p(1-ξ))^n) - ξ((1+p)^n-1))

from which the correlation may be found, by an appropriate scaling
using the variances given earlier.



  §.§.§ Binomial sampling

 
The generating function for Binomial sampling with n trials and
success probability q is


    (1-q + qt)^n



Let F(t,s) denote the joint for main T and stutter S peaks
for k cycles based on an initial single allele. Let N denote
the mean number of main alleles S the mean number of stutter
allele, and similarly for the variance N and S. Then for an
initial set of n alleles sampled with probability p the joint is given by


    Q(t,s)=(1-q+qF(t,s))^n


With abuse of notation, let t denote the mean number of main
alleles for this , , so that
t = ∂ Q/∂ t |_t=s=1, 

Thus
 

    ∂ Q(t,s)/∂
        t=nq(1-q+qF(t,s))^n-1∂ F(t,s)/∂ t

 
from which it follows that

    t = nq  N

 
Similarly
 
 
    ∂ Q(t,s)/∂
         s=nq(1-q+qF(t,s))^n-1∂ F(t,s)/∂ s

 
 from which it follows that
 
    s = nq  S

 
 Taking the second derivatives, we have
 
 
    ∂^2Q(t,s)/∂ t^2   =n(n-1)q^2(1-q+qF(t,s))^n-2(∂ F(t,s)/∂ t)^2 + nq(1-q+qF(t,s))^n-1∂^2 F(t,s)/∂ t^2
     t(t-1)    = n(n-1)q^2( N)^2 + nq ( N(N-1)
     t    =  n(n-1)q^2( N)^2 + nq ( N(N-1) + (nq  N) - (nq  N)^2
       = nq( N + (1-q) ( N)^2


 and similarly
 
    s = nq( S + (1-q) ( S)^2)


 Finally


    ∂^2Q(t,s)/∂ t∂ s   =nq(1-q+qF(t,s))^n-1∂ F(t,s)/∂ t∂ F(t,s)/∂ s
                                                         +
                                                         n(n-1)q^2(1-q+qF(t,s))^n-2∂^2F(t,s)/∂ t∂ s
     ts    = n(n-1)q^2 ( N  S) + nq ( NS)
    (t,s)    =  n(n-1)q^2 ( N  S) + nq ( NS) - (nq  N)(nq  S)
       =nq(  NS  - q( N)( S))
       = nq((N,S) + (1-q)( N)( S))


Hence


    (t,s)    = nq((N,S) + (1-q)( N)( S))/√(nq( N + (1-q) ( N)^2)√(nq( S + (1-q) ( S)^2)
       = (N,S) + (1-q)( N)( S))/√(( N + (1-q) ( N)^2)√(( S + (1-q) ( S)^2)




  §.§.§ Poisson sampling


 
Suppose that the number of alleles to be sampled has a Poison
distribution of mean λ. The using the for a Poisson
distribution, exp(λ(t-1)), we proceed in a similar manner as
for the Binomial sampling, using the joint 


    Q(t,s)=exp(λ(F(t,s)-1))



    ∂ Q(t,s)/∂ t   = λ∂ F(t,s)/∂ t Q(t,s)
     t    = λ N
    ∂ Q(t,s)/∂ s   = λ∂ F(t,s)/∂ s Q(t,s)
     s    = λ S
    
      ∂^2 Q(t,s)/∂ t^2   = λ^2 
                                                 (∂ F(t,s)/∂ t)^2 Q(t,s)
                                                 +
                                                 λ∂^2 F(t,s)/∂ t^2 Q(t,s)
    
      E t(t-1)    = λ^2 ( N)^2 + λ ( N(N-1))
    
      ∂^2 Q(t,s)/∂ s^2   = λ^2 
                                                 (∂ F(t,s)/∂ s)^2 Q(t,s)
                                                 +
                                                 λ∂^2 F(t,s)/∂ s^2 Q(t,s)
    
      E s(s-1)    = λ^2 ( S)^2 + λ ( S(S-1))
    
      ∂^2Q(t,s)/∂ t∂ s   =
                                                          λ^2 
                                                          ∂ F(t,s)/∂ t∂ F(t,s)/∂ sQ(t,s)
                                                          +
                                                          λ∂^2F(t,s)/∂ t∂ s Q(t,s)
     ts    = λ ( NS) + λ^2 ( N)( S)


from which it follows that

    (t,s)    = λ ( NS) + λ^2 ( N)( S) - λ^2 ( N)( S)
       = λ ( NS) 
     t    = λ^2 ( N)^2 + λ ( N(N-1)) + λ( N) - (λ N)^2
       = λ ( N^2)
     s    = λ ( S^2) ,
    (t,s)    = λ ( NS)/√(λ ( N^2))√(λ ( S^2))
       =  NS/√(( N^2)( S^2))

which does not depend on λ. In fact this result is the
limiting result of taking q→ 0 in the earlier Binomial sampling.




§ DERIVATION OF MOMENTS FOR GENOMIC MODEL, NO STUTTERS





 §.§ Moments of tagged amplicons

For convenience let t = t_a_d, and let D denote the derivative
with respect to t.  The marginal for the tagged amplicons is
given by the sum of such amplicons arising from the g and g_d
strand: we obtain this by multiplying the s.

    F(t) = G_n(1,1,1,t)G_d;n(1,1, t, 1),

from which

    DF(t) = [DG_n(1,1,1,t)]G_d;n(1,1, t, 1) + G_n(1,1,1,t)[DG_d;n(1,1, t, 1)]

hence the mean number of tagged amplicons is given by the sum of two
terms:

    E N_a_d = [DG_n(1,1,1,t)]_t=1 + [DG_d;n(1,1, t, 1)]_t=1

where we use G(1,1,1,) = G_d(1,1,1,1) = 1.  We may find each term on
the right hand side separately. Using the recurrence relation for the
joint of each genomic type, we may derive recurrence relations
for each of these two terms, which turn out to be coupled linear
equations that can be solved either by hand or computer algebra. We
consider each separately:



  §.§.§ DG

We may substitute 1 for every component, except t_a_d which we set
to t, in the joint recurrence relations. Assume this is done.
Then differentiating with respect to t and substituting t=1 we
obtain the following recurrence relations:
  

    DG_n+1    = (1-p_g)DG_n + p_g(DG_n + DH_d;n)
       = DG_n  + p_gDH_d;n
    
      DH_d;n+1    = (1-p_h_d)DH_d;n + p_h_d(DH_d;n + DA_n)
       = DH_d;n  + p_h_d DA_n
    
      DA_n+1    = (1-p_a)DA_n + p_a(DA_n +DA_d;n)
       = DA_n + p_a DA_d;n
    
      DA_d;n+1    = (1-p_a_d)DA_d;n + p_a_d(DA_d;n+DA_n)
       = DA_d;n  + p_a_dDA_n

with initial conditions DG_0 = DH_d;0 = DA_0 = 0 and
DA_d;0 = 1.

Now consider the last two coupled equations:

    DA_n+1    =DA_n + p_a DA_d;n
    
      DA_d;n+1    = DA_d;n  + p_a_dDA_n

We may write this in matrix form


    [   DA_n+1; DA_d;n+1 ]  =   [     1   p_a; p_a_d     1 ]  [   DA_n; DA_d;n ]

from which we clearly obtain



    [   DA_n; DA_d;n ]  =   [     1   p_a; p_a_d     1 ]^n
      [ 0; 1 ]

and hence

    DA_n =
      [ 1 0 ][     1   p_a; p_a_d     1 ]^n
      [ 0; 1 ].

Let us denote the square matrix (not raised to the power n) by P:

    P = [     1   p_a; p_a_d     1 ].

Then

    DA_n =
      [ 1 0 ]
      P^n
      [ 0; 1 ]

and from the earlier recurrence for DH_d;

    DH_d;n+1 = DH_d;n  + p_h_d DA_n

we obtain

    DH_d;n = p_h_d∑_j=0^n-1[ 1 0 ]
      P^j
      [ 0; 1 ]
      = p_h_d[ 1 0 ]P^n - I/P-I[ 0; 1 ]

From the earlier recurrence for DG;

    DG_n+1 = DG_n  + p_g DH_d;n

we therefore obtain

    DG_n = p_gp_h_d∑_j=0^n-1[ 1 0 ]P^j - I/P-I[ 0; 1 ]
      = p_gp_h_d[ 1 0 ]P^n -I - n(P-I)/(P-I)^2[ 0; 1 ]

which is our final result for the expected number of amplicons
arising from the g genomic strand.

Note that if all the amplification probabilities are equal to 1, we
obtain:

    P^n = [ 1 1; 1 1 ]^n =
      [ 2^n-1 2^n-1; 2^n-1 2^n-1 ]

which leads to DG_n = 2^n-1 - n obtained earlier.

We now need to do a similar calculation for the g_d strand.



  §.§.§ DG_d

We may substitute 1 for every component, except t_a_d which we set
to t, in the joint recurrence relations. Assume this is done.
Then differentiating with respect to t and substituting t=1 we
obtain the following recurrence relations:
  

    DG_d;n+1    = (1-p_g_d)DG_d;n + p_g_d(DG_d;n + DH_n)
       = DG_d;n  + p_g_dDH_n
    
      DH_n+1    = (1-p_h)DH_n + p_h(DH_n + DA_d;n)
       = DH_n  + p_h DA_d;n
    
      DA_n+1    = (1-p_a)DA_n + p_a(DA_n +DA_d;n)
       = DA_n + p_a DA_d;n
    
      DA_d;n+1    = (1-p_a_d)DA_d;n + p_a_d(DA_d;n+DA_n)
       = DA_d;n  + p_a_dDA_n

with initial conditions DG_0 = DH_d;0 = DA_0 = 0 and
DA_d;0 = 1. Note that the last two coupled equations are as before
for DG, and so have the same solution as given above:



    [   DA_n; DA_d;n ]  =   [     1   p_a; p_a_d     1 ]^n
      [ 0; 1 ]   =P^n
      [ 0; 1 ]

The equation for H reads:

    H_n+1= DH_n  + p_h DA_d;n

but

    DA_d;n =[ 0 1 ]
      P^n
      [ 0; 1 ]

Hence

    DH_n = p_h∑_j=0^n-1[ 0 1 ]
      P^j
      [ 0; 1 ]
      = p_h[ 0 1 ]P^n - I/P-I[ 0; 1 ]

Finally from the recurrence for DG_d:


    DG_d;n+1 = DG_d;n  + p_g_d DH_n

we obtain


    DG_d;n = p_g_dp_h∑_j=0^n-1[ 0 1 ]P^j - I/P-I[ 0; 1 ]
      = p_g_dp_h[ 0 1 ]P^n -I - n(P-I)/(P-I)^2[ 0; 1 ]

which is our final result for the expected number of amplicons
arising from the g_d genomic strand.

Note that is all amplification efficiencies are equal to 1, then we
obtain DG_d;n = 2^n-1-1 obtained earlier.



  §.§.§ DG and DG_d combined


The mean number of tagged amplicons is given by the sum of the two
results obtained above, which may be written in the form:


    [ p_gp_h_d,  p_g_dp_h ]P^n -I - n(P-I)/(P-I)^2[ 0; 1 ]



If p_gp_h_d= p_g_dp_h, which will happen if p_g= p_g_d
and p_h_d= p_h, this simplifies to:


    p_gp_h[ 1 1 ]P^n -I - n(P-I)/(P-I)^2[ 0; 1 ]

with

    P = [     1   p_a; p_a_d     1 ]

If p_a = p_a_d = p then

    P = [ 1 p; p 1 ]

then

    P^n = 1/2[   (1+p)^n+(1-p)^n (1+p)^n - (1-p)^n; (1+p)^n - (1-p)^n (1+p)^n + (1-p)^n ]

and

    (P-I)^-2 = 
    [ p^-2    0;    0 p^-2 ]
    = 1/p^2I

 
 After a little matrix algebra
 
 we obtain that the mean number of tagged amplicons is given by


    E A_d = p_gp_h/p^2( (1+p)^n - np -1)

 
Note that if p_g = p_h = p = 1 this reduces to 2^n-n-1 obtained
earlier for the Eulerian numbers.  If we have instead p_g = p_h = p
we obtain


    E A_d = (1+p)^n - np -1,


which can be compared to the value (1+p)^n of the simple model
that starts with a single amplicon.

Note that, if starting with m amplicons these mean values are simply
multiplied by m.

Returning to the more general result:


    [ p_gp_h_d,  p_g_dp_h ]P^n -I - n(P-I)/(P-I)^2[ 0; 1 ]

with

    P = [     1   p_a; p_a_d     1 ]

The eigenvalues of P are 1 + √(p_ap_a_d) and
1 - √(p_ap_a_d), and it can be shown that:


    P^n = 1/2[                     (1+√(p_ap_a_d))^n+(1-√(p_ap_a_d))^n                                                       ,     √(p_a/p_a_d)((1+√(p_ap_a_d))^n - (1-√(p_ap_a_d))^n); √(p_a_d/p_a)((1+√(p_ap_a_d))^n -
        (1-√(p_ap_a_d))^n)                                                       ,               (1+√(p_ap_a_d))^n +
        (1-√(p_ap_a_d))^n ]

and

    (P-I)^-2 =1/p_ap_a_d I

 
 More matrix algebra
 
 shows that


    (P^n -I - n(P-I))
      [ 0; 1 ]
      =
      [ √(p_a/p_a_d)[(1+√(p_ap_a_d))^n-(1-√(p_ap_a_d))^n])/2-np_a;                 [(1+√(p_ap_a_d))^n+(1-√(p_ap_a_d))^n]/2-1 ]


Hence the mean number of amplicons for the general case is given by
the row and column matrix inner product:



    [ p_gp_h_d/p_ap_a_d,  p_g_dp_h/p_ap_a_d ][ √(p_a/p_a_d)[(1+√(p_ap_a_d))^n-(1-√(p_ap_a_d))^n]/2-np_a;                [(1+√(p_ap_a_d))^n+(1-√(p_ap_a_d))^n]/2-1 ]

that is:

    p_gp_h_d/p_ap_a_d(
      √(p_a/p_a_d)(1+√(p_ap_a_d))^n-(1-√(p_ap_a_d))^n/2-np_a
    ) + p_g_dp_h/p_ap_a_d(
      (1+√(p_ap_a_d))^n+(1-√(p_ap_a_d))^n/2-1 )


If p_a = p_a_d = p, this simplifies to:


    p_gp_h_d/p^2( (1+p)^n-(1-p)^n/2-np ) +
    p_g_dp_h/p^2( (1+p)^n+(1-p)^n/2-1 )






  §.§.§ Variance of the number of amplicons


Recall that if F(t) is a for a random variable X then


    X    = dF(t)/dt|_t=1 = F'(1)
     X(X-1)    = dF(t)/dt|_t=1 = F”(1)
     X    = F”(1) + F'(1) - (F'(1))^2


We can apply this to

    F(t) = G_n(1,1,1,t)G_d;n(1,1, t, 1),

however is is simpler to consider the variances of the number of
tagged amplicons arising separately from the individual g and g_d
strands and add them: because they amplify independently their
variances will add to give the total variance of interest. This we now
do.

Recall that our basic vectorial has the form
  

    G_n+1   → G_n(1-p_g) + p_gG_nH_d;n
    
      H_d;n+1   → H_d;n(1-p_h_d) + p_h_dH_d;nA_n
    
      A_n+1   → A_n(1-p_a) + p_aA_nA_d;n
    
      A_d;n+1   → A_d;n(1-p_a_d) + p_a_dA_d;nA_n

with

    G_0    = t_g
    
      H_d,0   = t_h_d
    
      A_d   = t_a
    
      A_d,0   = t_a_d


Differentiating each twice with respect to t_a_d and setting all
the t's equal to 1 gives this set of linear recurrence relations


  

    D^2G_n+1   = D^2 G_n + 2p_gDG_nDH_d;n + p_g D^2H_d;n
    
      D^2H_d;n+1    =D^2H_d;n +2 p_h_dDH_d;nDA_n + p_h_dD^2A_n
    
      D^2A_n+1    = D^2A_n + 2p_aDA_nDA_d;n + p_aD^2A_d;n
    
      D^2A_d;n+1    = D^2 A_d;n +2 p_a_dDA_d;nDA_n + p_a_dD^2A_n


As before, we consider the last two coupled set of equations. These
may be written in matrix form as follows:



    [   D^2A_n+1; D^2A_d;n+1 ]  =   [     1   p_a; p_a_d     1 ][   D^2A_n; D^2A_d;n ]   +    2[   p_aDA_nDA_d;n; p_a_dDA_d;nDA_n ]

Previously we had the first order solution

    [   DA_n; DA_d;n ]  =   [     1   p_a; p_a_d     1 ]^n
      [ 0; 1 ]   =P^n
      [ 0; 1 ]

where

    P^n = 1/2[                     (1+√(p_ap_a_d))^n+(1-√(p_ap_a_d))^n                                                       ,     √(p_a/p_a_d)((1+√(p_ap_a_d))^n - (1-√(p_ap_a_d))^n); √(p_a_d/p_a)((1+√(p_ap_a_d))^n -
        (1-√(p_ap_a_d))^n)                                                       ,               (1+√(p_ap_a_d))^n +
        (1-√(p_ap_a_d))^n ]

hence

    [   DA_n; DA_d;n ]  =   1/2[ √(p_a/p_a_d)((1+√(p_ap_a_d))^n - (1-√(p_ap_a_d))^n);               (1+√(p_ap_a_d))^n + (1-√(p_ap_a_d))^n ]

and it follows that

    DA_nDA_d;n = 1/2√(p_a_d/p_a)(
      (1+√(p_ap_a_d))^2n - (1-√(p_ap_a_d))^2n)

This could be substituted into (<ref>) and the resulting
in-homogeneous difference equation may be solved (initial conditions
are that the second derivatives D^2A_0 = D^2A_d;0 = 0). It gets a
bit complicated, so it is simpler to work from the following scalar
equation and the substitute matrices for scalar.

So consider the recurrence relation

    y_n+1 = p y_n + c(a^n - b^n)

where y_0 = 0. Iteration gives the solution sequence:

    y_1    = 0
    
      y_2    = (a-b)c
    
      y_3    = p(a-b)c + (a^2-b^2)c
    
      y_4    = p^2(a-b)c + p(a^2-b^2)c + (a^3-b^3)c
       ⋮
    
      y_n    = ∑_i=0^n-2 p^n-i(a^i+1 - b^i+1)c
       = (ap^n-1-a^n-1/p-a - bp^n-1-b^n-1/p-b)c
  


Referring back , we may substitute
c = √(p_a_d/p_a)p_ap_a_d,
a = (1+√(p_a_dp_a))^2 I, b = (1-√(p_a_dp_a))^2I, where
I is the 2× 2 identity matrix, and p=P, to recover a matrix
formula for the solution of y_n ≡D^2A_nD^2A_d;n.

The result is a fairly complex matrix expression that has to be picked
apart to substitute various terms into the equations to solve for the
h and h_d, which themselves in turn have to be substituted into
equations to solve for g and g_d. This could possibly be done
using a computer algebra system, the results can be expected to be
messy - even more so when we introduce stutter - so we stop the
algebraic analysis here.

It seems appropriate therefore to generate matrix expressions at the
outset which would be more amenable to numerical evaluation. We now
follow this more direct approach. We could proceed by considering the
matrix formulation for g and g_d separately, or by using a
combined approach. We do the latter, they lead to the same results.




 §.§ Moments from matrix analysis



  §.§.§ First moment: Mean as matrix expression


The recurrence relation for the mean may be expressed as



    [   DG_n+1; DG_d;n+1;   DH_n+1; DH_d;n+1;   DA_n+1; DA_d;n+1 ]  =   [     1     0     0   p_g     0     0;     0     1 p_g_d     0     0     0;     0     0     1     0     0   p_h;     0     0     0     1 p_h_d     0;     0     0     0     0     1   p_a;     0     0     0     0 p_a_d     1;       ]  [   DG_n; DG_d;n;   DH_n; DH_d;n;   DA_n; DA_d;n ]   = P [   DG_n; DG_d;n;   DH_n; DH_d;n;   DA_n; DA_d;n ]


where P is now the 6 × 6 matrix. If we denote the column
matrix on the right by Y_n, then we have

    Y_n+1 = PY_n →  Y_n = P^n Y_0

where Y_0 is the transpose of the row vector (0,0,0,0,0,1), and
the mean number of tagged amplicons is given by the sum of the
elements in the first two rows of Y_n.

Numerically this is straightforward to extract given the values of the
various amplification probabilities.



  §.§.§ Second moment: as matrix expression


The matrix equation for the second moment follows the pattern of the
first order moment matrix equation.  Let Z_n denote the column
vector of second derivatives, that is


    Z_n =
      [   D^2G_n; D^2G_d;n;   D^2H_n; D^2H_d;n;   D^2A_n; D^2A_d;n ]


Then the recurrence equations for the second moments may be written as


    Z_n+1 = P Z_n + 2
      [   p_g DG_nDH_d;n; p_g_d DG_d;nDH_n;   p_h DH_nDA_d;n; p_h_d DH_d;nDA_n;   p_a DA_nDA_d;n; p_a_d DA_d;nDA_n ]
      = PZ_n + 2
      [    p_gY_n(1)Y_n(4); p_g_d Y_n(2)Y_n(3);   p_h Y_n(3)Y_n(6);  p_h_dY_n(4)Y_n(5);   p_a Y_n(5)Y_n(6); p_a_d Y_n(6)Y_n(5) ]


where Y_n(i) is the element in the i^th row of the column vector
Y_n. We have the initial condition Z_0 = (0,0,0,0,0,0)^T.

Thus we can iterate to solve sequentially for the pairs of column
vectors Y_1, Z_1; Y_2, Z_2; …; Y_n, Z_n.

Now Y_n(1) will denote the mean number of tagged amplicons arising
from the g strand, and Z_n(1) will denote the corresponding
expectation E M(M-1) for the number of tagged amplicons arising for
the g strand. Hence the variance of tagged amplicons arising form
the g strand will be

    Z_n(1) + Y_n(1) - Y^2_n(1)

Similarly the variance of tagged amplicons from the g_d strand will
be

    Z_n(2) + Y_n(2) - Y^2_n(2)

Hence the total variance of tagged amplicons will be

    Z_n(1) + Y_n(1) - Y^2_n(1) + Z_n(2) + Y_n(2) - Y^2_n(2)

If there are m genomic strands to begin with, we simply multiply
this variance expression by m to get the required total variance.




§ DERIVATION OF MOMENTS FOR GENOMIC MODEL, SINGLE STUTTERS






 §.§ Moments


The combined set of iterative equations of the branching process is:


    t_g    → t_g(1-p_g) + p_g(1-ξ)t_gt_h_d + p_gξ t_gt_h_sd
    
      t_g_d   → t_g_d(1-p_g_d) + p_g_d(1-ξ)t_g_dt_h + p_g_dξ t_g_dt_h_s
    
      t_h_d   → t_h_d(1-p_h_d) + p_h_d(1-ξ)t_h_dt_a+ p_h_dξ t_h_dt_a_s
    
      t_h_sd   → t_h_sd(1-p_h_sd) + p_h_sd t_h_sd  t_a_s
    
      t_h   → t_h(1-p_h) + p_h(1-ξ)t_ht_a_d+ p_hξ t_ht_a_sd
    
      t_h_s   → t_h_s(1-p_h_s) + p_h_s t_h_s t_a_s
    
      t_a   → t_a(1-p_a) + p_a(1-ξ)t_at_a_d + p_aξ t_at_a_sd
    
      t_a_d   → t_a_d(1-p_a_d) + p_a_d(1-ξ) t_a_dt_a + p_a_dξ t_a_dt_a_s
    
      t_a_s   → t_a_s(1-p_a_s) + p_a_st_a_st_a_sd
    
      t_a_sd   → t_a_sd(1-p_a_sd) + p_a_sdt_a_sdt_a_s

which lead directly to the coupled set of equations for the vectorial
generating function:

    G_n+1    = G_n[(1-p_g)+ p_g(1-ξ) H_n;d + p_gξ  H_n;sd]
    
      G_n+1;d   = G_n;d [ (1-p_g_d) + p_g_d(1-ξ)H_n + p_g_dξ H_n;s] 
    
      H_n+1    = H_n[ (1-p_h) + p_h(1-ξ)A_n;d+ p_hξ A_n;sd]
    
      H_n+1;d    = H_n;d[ (1-p_h_d) + p_h_d(1-ξ)A_n+ p_h_dξ A_n;s]
    
      H_n+1;s   =  H_n+1;s [(1-p_h_s) + p_h_sA_n;sd] 
    
      H_n+1;sd    = H_n;sd[(1-p_h_sd) + p_h_sd A_n;s] 
    
      A_n+1    = A_n[ (1-p_a) + p_a(1-ξ)A_n;d + p_aξ A_n;sd] 
    
      A_n+1;d    = A_n;d [ (1-p_a_d) + p_a_d(1-ξ) A_n + p_a_dξ A_n;s] 
    
      A_n+1;s    = A_n;s[(1-p_a_s) + p_a_sA_n;sd ]
    
      A_n+1;sd    = A_n;sd[(1-p_a_sd) + p_a_sdA_n;s]


Differentiating once for the mean, we obtain the matrix equation

    [    DG_n+1;  DG_n+1;d;    DH_n+1;  DH_n+1;d;  DH_n+1;s; DH_n+1;sd;    DA_n+1;  DA_n+1;d;  DA_n+1;s; DA_n+1;sd ]
      =
      [          1          0          0   p_g(1-ξ)          0       p_gξ          0          0          0          0;          0          1 p_g_d(1-ξ)          0     p_g_dξ          0          0          0          0          0;          0          0          1          0          0          0          0   p_h(1-ξ)          0       p_hξ;          0          0          0          1          0          0 p_h_d(1-ξ)          0     p_h_dξ          0;          0          0          0          0          1          0          0          0          0      p_h_s;          0          0          0          0          0          1          0          0     p_h_sd          0;          0          0          0          0          0          0          1   p_a(1-ξ)          0       p_aξ;          0          0          0          0          0          0 p_a_d(1-ξ)          1     p_a_dξ          0;          0          0          0          0          0          0          0          0          1      p_a_s;          0          0          0          0          0          0          0          0     p_a_sd          1;            ]
      [    DG_n;  DG_n;d;    DH_n;  DH_n;d;  DH_n;s; DH_n;sd;    DA_n;  DA_n;d;  DA_n;s; DA_n;sd ]



Denoting the 10 × 10 square matrix by P, and the column
matrix on the right by Y_n, then we have

    Y_n+1 = PY_n →  Y_n = P^n Y_0

where Y_0 is the transpose of the row vector
(0,0,0,0,0,0,0,0,0,1), and the mean number of tagged amplicons is
given by the sum of the elements in the first two rows of Y_n. This
will be equal to P^n[1,10] + P^n[2,10].







Let Z denote the column of second derivatives, then we have the
following in-homogeneous matrix recurrence relation that is readily
solved numerically by iteration:


    Z_n+1 = P Z_n + 2
      [   p_gDG_n[(1-ξ) DH_n:d + ξ DH_n;sd];  p_g_dDG_n;d[(1-ξ) DH_n + ξ DH_n;s];   p_hDH_n[ (1-ξ)DA_n;d + ξ DA_n;sd];  p_h_dDH_n;d[ (1-ξ)DA_n + ξ DA_n;s];                  p_h_sDH_n;sDA_n;sd;                 p_h_sdDH_n;sdDA_n;s;  p_aDA_n[ (1-ξ)D A_n;d + ξ DA_n;sd]; p_a_dDA_n;d[ (1-ξ) DA_n + ξ DA_n;s];                 p_a_sDA_n;s DA_n;sd;                p_a_sdDA_n;sd DA_n;s ]



and similarly to before, the total variance of tagged amplicons in
stutter position is given by

    Z_n(1) + Y_n(1) - Y^2_n(1) + Z_n(2) + Y_n(2) - Y^2_n(2)




 §.§ Covariance with main peak


The vectorial generating function may be used to find the mean and
variance of the main peak, and also the covariance between main peak
and stutter peak.  For the main peak, the mean height is given by
multiplying P^n into the column matrix that has zero everywhere
except for a 1 on the seventh row, and then adding the values in the
first two rows of the resulting column matrix.  The variance is
similarly found.

The second derivatives obey the following recurrence relation.

    [    DDsG_n+1;  DDsG_n+1;d;    DDsH_n+1;  DDsH_n+1;d;  DDsH_n+1;s; DDsH_n+1;sd;    DDsA_n+1;  DDsA_n+1;d;  DDsA_n+1;s; DDsA_n+1;sd ]
        = P
        [    DDsG_n;  DDsG_n;d;    DDsH_n;  DDsH_n;d;  DDsH_n;s; DDsH_n;sd;    DDsA_n;  DDsA_n;d;  DDsA_n;s; DDsA_n;sd ]
        +
        [   p_gDG_n[(1-ξ) DsH_n:d + ξ DsH_n;sd];  p_g_dDG_n;d[(1-ξ) DsH_n + ξ DsH_n;s];   p_hDH_n[ (1-ξ)DsA_n;d + ξ DsA_n;sd];  p_h_dDH_n;d[ (1-ξ)DsA_n + ξ DsA_n;s];                   p_h_sDH_n;sDsA_n;sd;                  p_h_sdDH_n;sdDsA_n;s;  p_aDA_n[ (1-ξ)Ds A_n;d + ξ DsA_n;sd]; p_a_dDA_n;d[ (1-ξ) DsA_n + ξ DsA_n;s];                  p_a_sDA_n;s DsA_n;sd;                 p_a_sdDA_n;sd DsA_n;s ]
        +
        [   p_gDsG_n[(1-ξ) DH_n:d + ξ DH_n;sd];  p_g_dDsG_n;d[(1-ξ) DH_n + ξ DH_n;s];   p_hDsH_n[ (1-ξ)DA_n;d + ξ DA_n;sd];  p_h_dDsH_n;d[ (1-ξ)DA_n + ξ DA_n;s];                  p_h_sDsH_n;sDA_n;sd;                 p_h_sdDsH_n;sdDA_n;s;  p_aDsA_n[ (1-ξ)D A_n;d + ξ DA_n;sd]; p_a_dDsA_n;d[ (1-ξ) DA_n + ξ DA_n;s];                 p_a_sDsA_n;s DA_n;sd;                p_a_sdDsA_n;sd DA_n;s ]





§ A FURTHER LOOK AT TARGET AND STUTTER PEAK HEIGHT CORRELATIONS






In <cit.> the authors looked at developing a
model for the stutter ratio, defined as the ratio of the stutter peak
to main peak, such that the mean stutter ratio is linearly dependent
on the longest uninterrupted sequence (LUS) of repeats. They used
controlled experimental single-source samples, 289 in all, with a
target of 1ng of DNA and amplified with a 25 RFU detection
limit. After discarding loci which had heterozygous loci one repeat
apart, 2323 heterozygous loci were left for analysis.


It is important to note that they define stutter ratio as the ratio of
the stutter peak height to the main peak height:

    SR_a = O_a-1/O_a

They define the total allelic product to be

    T_a = O_a-1+ O_a


They propose a lognormal distribution for the ratio of observed to
expected peak height for main and stutter peaks (they comment that a
gamma model is not suitable because of the heavy tails). From their
data they conclude that there is little evidence of correlation
between the main and stutter peaks given the target amount of
DNA. They write: “Unexpectedly the scatter plots
  in Fig. 8a and b indicate that there is no detectable correlation
  between stutter and allele in this biological model”. They found a
low correlation of around 0.11.

Now clearly having little to no correlation is quite at variance with
the branching model.  We suggest that this is an artefact of
grouping the data together within and across loci. First we recast
their results. In their Appendix 2 they give the following table for
the linear regression of their stutter ratio model (using the NGM
SELect system), in which the stutter ratio of an allele i, having  LUS L_i,  from a locus is equal to
SR_i = a + b× L_i:





  
        Locus     Intercept a     Slope b 
 
    1     D10S1248     -0.0576     0.0089 
 
    2     D12S391     -0.0571     0.0107 
 
    3     D16S539     -0.0502     0.0088 
 
    4     D18S51     -0.0297     0.0066 
 
    5     D19S433     -0.0302     0.0074 
 
    6     D1S1656     -0.0699     0.0106 
 
    7     D21S11     -0.0079     0.0059 
 
    8     D22S1045     -0.0881     0.0139 
 
    9     D2S1338     -0.0073     0.0062 
 
    10     D2S441     0.0004     0.0031 
 
    11     D3S1358     -0.0455     0.0092 
 
    12     D8S1179     -0.0148     0.0062 
 
    13     FGA     -0.0344     0.0066 
 
    14     SE33     0.0129     0.0041 
 
    15     TH01     -0.0208     0.0052 
 
    16     vWA     -0.0354     0.0078 
 
    



Now from the simple amplicon model, with amplification probability p
and conditional stutter probability ξ, the expected number of
amplicons heights for main and stutter peaks given k cycles and
starting from n_0 amplicons, are given by


    ET    = n_0(1 +p(1-ξ))^k
    
      ES    = n_0[(1 +p)^k -  (1 +p(1-ξ))^k]


Thus we can approximate, with ξ dependent on L:


    SR = a + b L ≈ES/ET =  (1+p/(1+p(1-ξ)))^k - 1
    ≈kp/1+pξ + k(k-1)p^2/2(1+p)^2ξ^2 +
    0(ξ^3)


From this we obtain to a good approximation, by taking the linear term
only, the incremental change in mean stutter ratio by a unit increase
in LUS:

    Δξ = b1+p/kp


For p ∈ [0.8,1], (1+p)/p varies between 2 and
2.25, hence to a good approximation Δξ = 2b/k.


So now we can collect results. Given the nature of the experimental
setup, with a target DNA amount of 1ng, (corresponding to around 152
cells) we may take the number of alleles from a heterozygous marker to
be Pois(152). In addition the LUS in their plots varies from 10-15
for most markers, some having a larger range, some smaller.  So for
simplicity we shall simulate LUS values from the same range for each
of the 16 markers, (also corresponding therefore to a uniform
distribution of alleles).  The following table collects individual
marker correlations, based on around 2000 simulations on each marker,
under various scenarios, but using the ξ values for each marker
using the SR formula above; we take k=28 and p=0.85 in all
simulations.




  
    Locus    
           LUS=12   
                   LUS=12   
                           LUS∼Unif[8,15]   
                                                LUS∼Unif[8,15] 
 
            n_0=152   
                     n_0∼ Pois(152)   
                                          n_0=152   
                                                     n_0∼ Pois(152) 
 D10S1248    0.2005406    0.7255809    -0.5541554    -0.07269729 

    D12S391    0.2290437    0.773266    -0.6050955    -0.04926867 

    D16S539    0.2054748    0.7379974    -0.5366395    -0.03364616 

    D18S51    0.1665038    0.7108246    -0.4260814    0.06089547 

    D19S433    0.248041    0.7182632    -0.4844209    0.06972643 

    D1S1656    0.1848627    0.7401926    -0.6481533    -0.0976838 

    D21S11    0.1952987    0.7507546    -0.3596321    0.2097556 

    D22S1045    0.217241    0.7738013    -0.7170921    -0.1489542 

    D2S1338    0.2145392    0.7708068    -0.3787377    0.1656755 

    D2S441    0.1509088    0.6713596    -0.1463518    0.2585785 

    D3S1358    0.2333797    0.7477239    -0.5480802    -0.03041745 

    D8S1179    0.2451524    0.731677    -0.3884186    0.1381175 

    FGA    0.1308529    0.6779035    -0.4587745    0.06434562 

    SE33    0.1765904    0.6558634    -0.2304545    0.1803573 

    TH01    0.1666116    0.6661396    -0.3677087    0.1420518 

    vWA    0.1850046    0.7455873    -0.5083684    0.0339069 
 



We see that keeping LUS and n_0 fixed we get moderate correlations
as in earlier chapter. Making the starting number of alleles Poisson
distributed and keeping LUS fixed leads to high correlation between
the stutter and peak amplicon numbers, also as found earlier. Now
making LUS random (to simulate a range of genotypes) but keeping n_0
fixed the correlations are negative. Combining the randomization of
LUS and making n_0∼ Pois(132) yields the correlations in the
final column, which are largely quite small with some positive and
some negative.


Finally we do a simulation in which we aggregate the various loci
simulations: here are 16 runs of such a simulation which mirrors the
experimental data of <cit.>. We see that all the
correlations are negative, but very close to zero.



   -0.003916072    -0.02691701    -0.0257117    -0.03113958 

    -0.02115133    -0.02109356    -0.02308784    -0.05802165 

    -0.02720094    -0.0288306    -0.0345473    -0.07301637 

    -0.01966967    -0.031646    -0.0649325    -0.02202189 
 



The following scatterplot was obtained from one simulation, in which
we plot the log of the sumulated number of amplicons divided by their
expected number (itself estimated from simulation) for the target and
stutter alleles, so that it corresponds to the quantities plotted in
Figure 8b of <cit.>. The plot below has a smaller
range than their Figure 8b, which could be because we have simulated
from the amplicon model rather than the genomic model.  The sample
correlation of the points in the plot is 0.157, close to the figure of
0.11 they found from their experimental data.




  
    < g r a p h i c s >




Hence, on the basis of this
simulation,  we deduce that the bivariate  branching process model showing high main allele 
and stutter-peak correlations on individual alleles is
compatible with the experimental results of 
<cit.> which aggregates observations of  alleles within and between loci.


Support for this may be found in a dataset of experimental single
source and mixed profiles released by Boston University in
2012[see
  <http://www.bu.edu/dnamixtures/pages/help/introduction/>] These
were samples prepared from extracted DNA from 4 anonymous individuals
(labelled A,B,C and D) diluted to prepare controlled DNA amounts, and
amplified for four kits. In the single source samples there many
independent amplifications of each individual's DNA, allowing the
target-stutter peak height data to be found for individual alleles for
a given amount of initial DNA, making sure that neither are
"contaminated" by neighbouring alleles.  Looking at the genotypes , we
see that for person A allele 25 on the locus D2D1338 is a possible
target allele. Taking all the single source samples on 1ng and 10
seconds injection time for person A amplified with the
Identifiler™kit (chosen here because it yields the most
number of data points), the author extracted target (25) and stutter
(24) peak height values to produce the following scatterplot , having
a sample correlation of 0.929.





    < g r a p h i c s >




Bearing in mind the way the samples were prepared, we should expect
Poisson sampling for the initial number of amplicons prior to ,
and hence expect a high correlation to result as found. Looking at the
stutter proportions (stutter/total) which average around 0.08 for the
11 data points, we see that the results are very consistent with branching modelling approach for stutter.

In conclusion, we suggest that the combined effects of aggregating
across loci and alleles might be giving a misleading impression that
there is no significant correlation.  Dedicated experiments
concentrating on specific alleles and their stutter product, that do
not aggregate data over different alleles either within or between
loci, could be carried out to confirm or refute this.



§ SAMPLE CODE FOR GENERATING DISTRIBUTIONS AND PLOTS



 §.§ Python scripts





  §.§.§ code for fig:ampkd.







  §.§.§ code for fig:ampfft28







  §.§.§ code for
  alg-ampliconjointfft
 


[Joint distribution for target and stutter amplicons:
    Python code]
  

  
Note that even for M=1 and with a typical K=28 cycles in , the
size of the 2-dimensional array F to hold the joint distribution
will be N=2^56≈ 7.2× 10^16, which is far in excess
of the memory of any computer.  To fill the array with double
precision values would require approximation 4.2 billion Gigabytes of
Ram. Filling the array would also take an excessive amount of
time. Hence this should be used only for small K and M values.



  §.§.§ code for
  alg-ampliconstuttermargfft
 

[Marginal distribution for stutter amplicons: Python
    code]
  

  



 §.§ R scripts




  §.§.§ code for generating target marginal distribution, no
  stutters
 

R code for the distribution of the total number of amplicons arising
from K amplification cycles, amplification probability p on each
cycle, starting with M amplicons binomially sampled prior to
amplification with probability ϕ.  Note that R vector indices
start at 1 instead of 0, hence on line 3 we set F[2] = 1 (rather
than F[1]=1).

[Amplicon probability distribution using R]
  







  §.§.§ code to generate
  fig:dropoutamp





  §.§.§ code for generating marginal distributions for the
  genomic model, implementing alg:gemomemarg

 





 §.§ Julia scripts




  §.§.§ code for generating fig:fulltraunc

Julia, like R, uses arrays with indices starting from 1.  You will
need to install the PyPlot package to run this script, and will
require approximately 10Gb of free ram to do the full evaluation
with K=28.





tocsection



50


urlstyle

[Alfonse et al.(2016)Alfonse, Garrett, Swaminathan, Peters, Wellner,
  Yearwood-Garcia, Taranow, Lun, Duffy, Médard,
  et al.]alfonse2016development
L. E. Alfonse, A. D. Garrett, H. Swaminathan, K. C. Peters, G. Wellner,
  X. Yearwood-Garcia, L. M. Taranow, S. Lun, K. R. Duffy, M. Médard, et al.
The development and release of a collection of computational tools
  and a large-scale empirical data set for validation: The PROVEDIt
  initiative.
In International Symposium on Human Identification, 2016.

[Alfonse et al.(2018)Alfonse, Garrett, Lun, Duffy, and
  Grgicak]alfonse2018large
L. E. Alfonse, A. D. Garrett, D. S. Lun, K. R. Duffy, and C. M. Grgicak.
A large-scale dataset of single and mixed-source short tandem repeat
  profiles to inform human identification strategies: PROVEDIt.
Forensic Science International: Genetics, 32:0 62–70,
  2018.

[Balding(2013)]balding:13
D. Balding.
Evaluation of mixed-source, low-template DNA profiles in forensic
  science.
Proceedings of the National Academy of Sciences of the United
  States of America, 1100 (30):0 12241–12246, 2013.

[Balding(2005)]balding:05
D. J. Balding.
Weight-of-evidence for Forensic DNA Profiles.
Statistics in Practice. Wiley, Chichester, UK, 2005.

[Balding and Steele(2015)]balding2015weight
D. J. Balding and C. D. Steele.
Weight-of-evidence for Forensic DNA Profiles.
John Wiley & Sons, 2015.

[Bleka et al.(2016)Bleka, Storvik, and Gill]bleka2016euroformix
Ø. Bleka, G. Storvik, and P. Gill.
Euroformix: an open source software based on a continuous model to
  evaluate str dna profiles from a mixture of contributors with artefacts.
Forensic Science International: Genetics, 21:0 35–44,
  2016.

[Bleka et al.(2017)Bleka, Eduardoff, Santos, Phillips, Parson, and
  Gill]bleka2017open
Ø. Bleka, M. Eduardoff, C. Santos, C. Phillips, W. Parson, and P. Gill.
Open source software euroformix can be used to analyse complex SNP
  mixtures.
Forensic Science International: Genetics, 31:0
  105–110, 2017.

[Bodner et al.(2016)Bodner, Bastisch, Butler, Fimmers, Gill,
  Gusmão, Morling, Phillips, Prinz, Schneider,
  et al.]bodner2016recommendations
M. Bodner, I. Bastisch, J. M. Butler, R. Fimmers, P. Gill, L. Gusmão,
  N. Morling, C. Phillips, M. Prinz, P. M. Schneider, et al.
Recommendations of the dna commission of the international society
  for forensic genetics (isfg) on quality control of autosomal short tandem
  repeat allele frequency databasing (strider).
Forensic Science International: Genetics, 24:0
  97–102, 2016.

[Bright et al.(2013)Bright, Taylor, Curran, and
  Buckleton]bright2013developing
J.-A. Bright, D. Taylor, J. M. Curran, and J. S. Buckleton.
Developing allelic and stutter peak height models for a continuous
  method of DNA interpretation.
Forensic Science International: Genetics, 70
  (2):0 296–304, 2013.

[Butler(2005)]book:Butler
J. M. Butler.
Forensic DNA Typing: Biology, Technology, and Genetics of STR
  Markers, Second Edition.
Elsevier Academic Press, 2005.
ISBN 0121479528.

[Butler(2011)]butler2011advancedmethod
J. M. Butler.
Advanced topics in forensic DNA typing: methodology.
Academic Press, 2011.

[Butler(2014)]butler2014advancedinterp
J. M. Butler.
Advanced topics in forensic DNA typing: interpretation.
Academic Press, 2014.

[Butler et al.(2003)Butler, Schoske, Vallone, Redman, and
  Kline]butler:etal:03
J. M. Butler, R. Schoske, P. M. Vallone, J. W. Redman, and M. C. Kline.
Allele frequencies for 15 autosomal loci on U.S. Caucasian,
  African American, and Hispanic populations.
Journal of Forensic Science, 480 (4):0
  908–911, 2003.
Available online at www.astm.org.

[Cowell(2016)]cowell2016combining
R. Cowell.
Combining allele frequency uncertainty and population substructure
  corrections in forensic dna calculations.
Forensic Science International: Genetics, 23:0
  210–216, 2016.

[Cowell et al.(2015)Cowell, Graversen, Lauritzen, and
  Mortera]cowell2015analysis
R. Cowell, T. Graversen, S. Lauritzen, and J. Mortera.
Analysis of forensic DNA mixtures with artefacts.
Journal of the Royal Statistical Society: Series C (Applied
  Statistics), 640 (1):0 1–48, 2015.

[Cowell(2009)]cowell:09
R. G. Cowell.
Validation of an STR peak area model.
Forensic Science International: Genetics, 30
  (3):0 193–199, 2009.
ISSN 1872-4973.
10.1016/j.fsigen.2009.01.006.
URL
  <http://www.sciencedirect.com/science/article/pii/S1872497309000088>.

[Duffy et al.(2016)Duffy, Gurram, C Peters, Wellner, and
  Grgicak]duffy:etal:2016
K. Duffy, N. Gurram, K. C Peters, G. Wellner, and C. Grgicak.
Exploring str signal in the single- and multi-copy number regimes:
  Deductions from an in silico model of the entire dna laboratory process.
Electrophoresis, 38, 12 2016.

[Egeland et al.(2003)Egeland, Dalen, and Mostad]egeland2003estimating
T. Egeland, I. Dalen, and P. F. Mostad.
Estimating the number of contributors to a DNA profile.
International journal of legal medicine, 1170
  (5):0 271–275, 2003.

[Egeland et al.(2015)Egeland, Kling, and
  Mostad]egeland2015relationship
T. Egeland, D. Kling, and P. Mostad.
Relationship inference with familias and R: statistical methods
  in forensic genetics.
Academic Press, 2015.

[Evett and Weir(1998)]evett1998interpreting
I. W. Evett and B. S. Weir.
Interpreting DNA evidence: statistical genetics for forensic
  science.
Sunderland, Massachusetts, 1998.

[Findlay et al.(1997)Findlay, Frazier, Taylor, and
  Urquhart]Findlay1997
I. Findlay, R. Frazier, A. Taylor, and A. Urquhart.
Single cell DNA fingerprinting for forensic applications.
Nature, 389:0 555–556, 1997.

[Gill et al.(2005)Gill, Curran, and Elliot]gill:etal:2005
P. Gill, J. Curran, and K. Elliot.
A graphical simulation model of the entire DNA process associated
  with the analysis of short tandem repeat loci.
Nucleic Acids Research, 330 (2):0 632–643,
  2005.

[Good(1955)]good1955joint
I. Good.
The joint distribution for the sizes of the generations in a cascade
  process.
In Mathematical Proceedings of the Cambridge Philosophical
  Society, volume 51, pages 240–242. Cambridge University Press, 1955.

[Good(1949)]good:1949
I. J. Good.
The number of individuals in a cascade process.
Proceedings of the Cambridge Philosophical Society,
  45:0 360–363, 1949.

[Good(1950)]good1950
I. J. Good.
Probability and the Weighing of Evidence.
C. Griffin London, 1950.

[Good(1979)]good:79
I. J. Good.
Studies in the History of Probability and Statistics. XXXVII
  A. M. Turing's statistical work in World War II.
Biometrika, 660 (2):0 393–396, 1979.

[Graversen and Lauritzen(2015)]graversen2015computational
T. Graversen and S. Lauritzen.
Computational aspects of DNA mixture analysis.
Statistics and Computing, 250 (3):0 527–541,
  2015.

[Grisedale and van Daal(2014)]grisedale2014method
K. Grisedale and A. van Daal.
Method summary.
Biotechniques, 560 (3):0 145–147, 2014.

[Haned(2013)]forensim
H. Haned.
forensim: Statistical tools for the interpretation of forensic
  DNA mixtures, 2013.
URL <https://CRAN.R-project.org/package=forensim>.
R package version 4.3.

[Haned et al.(2011)Haned, Pene, Lobry, Dufour, and
  Pontier]haned2011estimating
H. Haned, L. Pene, J. R. Lobry, A. B. Dufour, and D. Pontier.
Estimating the number of contributors to forensic DNA mixtures:
  does maximum likelihood perform better than maximum allele count?
Journal of forensic sciences, 560 (1):0
  23–28, 2011.

[Heideman et al.(1984)Heideman, Johnson, and Burrus]heideman1984gauss
M. Heideman, D. Johnson, and C. Burrus.
Gauss and the history of the fast fourier transform.
IEEE ASSP Magazine, 10 (4):0 14–21, 1984.

[Jeffreys et al.(1985)Jeffreys, Wilson, and Thein]jeffreys:85
A. J. Jeffreys, A. J. Wilson, and S. L. Thein.
Individual specific `fingerprints' of human DNA.
Nature, 316:0 76–79, 1985.

[Lauritzen and Mortera(2002)]lauritzen2002bounding
S. L. Lauritzen and J. Mortera.
Bounding the number of contributors to mixed DNA stains.
Forensic science international, 1300 (2):0
  125–126, 2002.

[Lindley(1977)]lindley1977
D. V. Lindley.
A problem in forensic science.
Biometrika, 640 (2):0 207–213, 1977.

[Perlin and Szabady(2001)]article:Perlindata
M. Perlin and B. Szabady.
Linear mixture analysis: a mathematical approach to resolving mixed
  DNA samples.
Journal of Forensic Science, 46:0 1372–1378, 2001.

[Puch-Solis(2014)]puch2014dropin
R. Puch-Solis.
A dropin peak height model.
Forensic Science International: Genetics, 11:0 80–84,
  2014.

[R Core Team(2016)]Rmanual
R Core Team.
R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing, Vienna, Austria, 2016.
URL <https://www.R-project.org/>.

[Rao et al.(2011)Rao, Kim, and Hwang]rao2011fast
K. R. Rao, D. N. Kim, and J. J. Hwang.
Fast Fourier Transform-Algorithms and Applications.
Springer Science & Business Media, 2011.

[Scientific(2012)]scientific2012ampflstr
T. F. Scientific.
Ampflstr® identifiler® plus pcr
  amplification kit user’s guide.
Waltham, MA, USA: Thermo Fisher Scientific, 2012.

[Stolovitzky and Cecchi(1996a)]article:StolovitzkyCecchi
G. Stolovitzky and G. Cecchi.
Efficiency of DNA replication in the polymerase chain reaction.
Proceedings of the National Academy of Sciences, 93:0
  12947–12952, 1996a.

[Stolovitzky and Cecchi(1996b)]stolovitzky1996efficiency
G. Stolovitzky and G. Cecchi.
Efficiency of DNA replication in the polymerase chain reaction.
Proceedings of the National Academy of Sciences, 930
  (23):0 12947–12952, 1996b.

[Sun(1995)]article:SunPCR
F. Sun.
The polymerase chain reaction and branching processes.
Journal of Computational Biology, 20 (1):0
  63–86, 1995.

[Swaminathan et al.(2015)Swaminathan, Grgicak, Medard, and
  Lun]swaminathan2015nocit
H. Swaminathan, C. M. Grgicak, M. Medard, and D. S. Lun.
Nocit: A computational method to infer the number of contributors to
  DNA samples analyzed by STR genotyping.
Forensic Science International: Genetics, 16:0
  172–180, 2015.

[Tvedebrink et al.(2009)Tvedebrink, Eriksen, Mogensen, and
  Morling]tvedebrink2009estimating
T. Tvedebrink, P. S. Eriksen, H. S. Mogensen, and N. Morling.
Estimating the probability of allelic drop-out of STR alleles in
  forensic genetics.
Forensic Science International: Genetics, 30
  (4):0 222–226, 2009.

[Tvedebrink et al.(2010)Tvedebrink, Eriksen, Mogensen, and
  Morling]Tvedebrink2010
T. Tvedebrink, P. S. Eriksen, H. S. Mogensen, and N. Morling.
Evaluating the weight of evidence by using quantitative short tandem
  repeat data in DNA mixtures.
Applied Statistics, 590 (5):0 855 – 874,
  2010.

[Tvedebrink et al.(2012a)Tvedebrink, Eriksen, Asplund,
  Mogensen, and Morling]tvedebrink2012allelic
T. Tvedebrink, P. S. Eriksen, M. Asplund, H. S. Mogensen, and N. Morling.
Allelic drop-out probabilities estimated by logistic
  regression—further considerations and practical implementation.
Forensic Science International: Genetics, 60
  (2):0 263–267, 2012a.

[Tvedebrink et al.(2012b)Tvedebrink, Eriksen, Mogensen,
  and Morling]tvedebrink:deconv:12
T. Tvedebrink, P. S. Eriksen, H. S. Mogensen, and N. Morling.
Identifying contributors of DNA mixtures by means of quantitative
  information of STR typing.
Journal of Computational Biology, 190 (7):0
  887–902, 2012b.
10.1089/cmb.2010.0055.

[Walsh et al.(1992)Walsh, Erlich, and Higuchi]walsh1992preferential
P. S. Walsh, H. A. Erlich, and R. Higuchi.
Preferential pcr amplification of alleles: mechanisms and solutions.
Genome Research, 10 (4):0 241–250, 1992.

[Wang et al.(2006)Wang, Xue, and Birdwell]wang2006least
T. Wang, N. Xue, and J. D. Birdwell.
Least-square deconvolution: A framework for interpreting short tandem
  repeat mixtures.
Journal of Forensic Sciences, 510 (6):0
  1284–1297, 2006.

[Weusten and Herbergs(2012)]weusten2012stochastic
J. Weusten and J. Herbergs.
A stochastic model of the processes in PCR based amplification of
  STR DNA in forensic applications.
Forensic Science International: Genetics, 60
  (1):0 17–25, 2012.
URL <doi:10.1016/j.fsigen.2011.01.003>.


 
 



 
 




