=Paper=
{{Paper
|id=Vol-2387/20190017
|storemode=property
|title=Simulation as a Tool to Identify Dynamical Typology of Water Frog Hemiсlonal Population Systems
|pdfUrl=https://ceur-ws.org/Vol-2387/20190017.pdf
|volume=Vol-2387
|authors=Dmytro Shabanov,Marina Vladymyrova,Anton Leonov,Olga Biriuk,Marina Kravchenko,Quentin Mair,Olena Meleshko,Julian Newman,Olena Usova,Grygoriy Zholtkevych
|dblpUrl=https://dblp.org/rec/conf/icteri/ShabanovVLBKMMN19
}}
==Simulation as a Tool to Identify Dynamical Typology of Water Frog Hemiсlonal Population Systems==
<pdf width="1500px">https://ceur-ws.org/Vol-2387/20190017.pdf</pdf>
<pre>
Simulation as a Tool to Identify Dynamical Typology of
     Water Frog Hemiсlonal Population Systems

  Dmytro Shabanov1, Marina Vladymyrova1, Anton Leonov2, Olga Biriuk1,
  Marina Kravchenko1, Quentin Mair3, Olena Meleshko4, Julian Newman3,
               Olena Usova1 and Grygoriy Zholtkevych1
                 1 V.N. Karazin Kharkiv National University, Ukraine
                                2 AppScale Systems Ink.
                   3 Glasgow Caledonian University, Scotland, UK
  4 University Museum, Norwegian University of Science and Technology, Norway

                              shabanov@karazin.ua



   Abstract. Some related species give rise to interspecies hybrids with hemiclonal
   inheritance. The gametes of such hybrids transfer the set of hereditary infor-
   mation of one of the parental species. The water frog, Pelophylax esculentus, is
   an example of such hybrids. The hemiclonal hybrids together with their parental
   species form a biosystem for which the suggested name is Hemiclonal Population
   System (HPS). The phenomenon of interspecific hemiclonal reproduction of wa-
   ter frogs has been intensively explored for several decades, but insufficient study
   has been devoted to the mechanisms of the composition constancy and ecological
   stability of their population systems.
       In this paper we focus on sustainability and possible results of transformations
   (in terms of different genetic forms and of population dynamics of each of these
   forms) of HPSs that consist of diploid representatives. By means of a simulation
   model we evaluated the long-term consequences of parameters that were derived
   from empirical research on a natural HPS at a particular brief time-period. This
   research was carried out in the region of eastern Ukraine, the so-called Siverskyi
   Donets center of diversity of the Pelophylax esculentus complex. Our study de-
   scribes the development of a dynamic typology of the HPS by principal compo-
   nent analysis of data generated by the simulations. We investigate the space of
   the possible states of the HPS. The simulations helped to split this state-space
   into six areas of stability, each of which corresponds to a different type of stabil-
   ity. Before conducting the simulation study we assumed there were only four
   truly stable states. Two new states were identified as a result of using this model.

   Keywords: Dynamic Typology, Hemiclonal Inheritance, Pelophylax esculentus
   complex, Simulation Modelling.
1      Introduction

1.1    Motivation
The processes of transformation of natural biological systems (populations, ecosystems
etc.) are an important challenge for modern science. In many cases these processes are
not readily available for direct study (e.g. because of their long duration). In such cases
we have to reconstruct their transformation models based on (1) the observed diversity
of their natural states, (2) our hypothesis of what controls the change of such states, and
(3) available empirical data. A simulation model plays a key role in such studies.
Thanks to it, we can determine the possible directions of transformation of the target
systems and build their typology. One of the insufficiently studied types of biological
systems are those that create species capable of interspecies hybridization and their
hemiclonal hybrids. Reproduction of hemiclonal hybrids occurs in biological systems
which have been described as Hemiclonal Population Systems (HPSs) [15]. Water frogs
are a typical model organism for studying such systems.
    This paper describes a computer simulation study of Hemiclonal Population Systems
(HPS) of water frogs. First we introduce the peculiar characteristics of these biosys-
tems. In the next section we indicate the advantages of computer simulation as a mod-
elling tool together with the specific characteristics of the model presented here. Com-
puter simulation has been described as the collaboration between experimenting and
modelling [7]. The existing empirical data was used to determine model parameters and
to assess the performance of the model against real-world situations. We describe ex-
periments conducted upon the model itself and the way in which these can supplement
fieldwork findings.

1.2    Reproduction within HPSs
Most life forms which arise from sexual reproduction exhibit the typical biosystem hi-
erarchy in which organisms exist as part of populations that form species. In these pop-
ulations, genetically unique individuals produce sex cells (gametes), which bear a
unique genome (holistic unitary complex of hereditary information) resulting from re-
combination of two parental genomes (Fig.1 A). The consequence of this is the exist-
ence of a population gene pool (common pool of genes). Populations of organisms with
clonal reproduction, on the other hand, consist of clones, i.e. sets of genetically identical
organisms (with accuracy limited by the error rate of copying during reproduction).
These populations are presented as a set of relatively isolated lines consisting of genet-
ically identical maternal and child individuals.
   This work is dedicated to a relatively rare type of population reproduction which
differs from the two mentioned options. It is characteristic of hybridogenic complexes
of species. This complex consists of two parental frog species [11]. Their full zoological
names are Pelophylax lessonae (Camerano, 1882) and Pelophylax ridibundus (Pallas,
1771). They cross to produce hybrids named analogous to the species' name, e.g. Pelo-
phylax esculentus (Linnaeus, 1758).
Fig. 1. A comparison of the life cycle with fertilization and meiosis, which is characteristic of
organisms with sexual reproduction (A), and modification of such a cycle, which is characteristic
of interspecific hybrids with hemiclonal inheritance (B).

The differences between the genomes of the parental species lead to the fact that pro-
duction of normal gametes through recombination becomes impossible. Breeding hy-
brids is made possible by specific changes in their life cycle (Fig. 1B). One parent ge-
nome of such hybrids is eliminated from the germline cells (i.e. those cells which sub-
sequently form the sex cells, gametes). Thereafter, the second parental genome under-
goes endoreduplication (i.e. doubles without cell division). This gives rise to cells with
two genomes that are identical (with accuracy limited by the errors of copying). They
form genetically identical gametes.
   Hemiclonal hybrids differ both from usual organisms with recombinant reproduction
and from clonal organisms. A HPS, where hemiclonal hybrids are reproduced, differs
from normal recombinant and from clonal populations. The features of a HPS are as
follows [15]:

— Cooperative reproduction of individuals that differ in species composition of their
 genomes (i.e. representatives of parental species and hybrids of various genome
 compositions);
— Vertical transmission of the clonal genome lines, which can combine with other
 recombinant or clonal genomes;
— In the HPS, which includes individuals of the parental species, these individuals
 support the existence of a pool of recombining genes (gene pool), which corresponds
 to the gene pool of conventional monospecific populations;
— Cases of limited interspecies recombination are observed (i.e. transferring of frag-
 ments of genetic information from the genome of one parental species into the gene
 pool of another parental species).

    The most common typology of a HPS is associated with an indication of the parental
species and diploid hybrids and/or triploid hybrids, which are included in their compo-
sition. The presence of P. esculentus is indicated by the letter E, letter L stands for
P. lessonae, and letter R for P. ridibundus. The presence of polyploid P. esculentus is
denoted by the letter p. Thus, L-E-HPS consists of P. lessonae and diploid P. esculen-
tus, and R-E-HPS consists of P. ridibundus and diploid P. esculentus. Clonality of a
genome is designated by putting its symbol in brackets. The water frogs' sex determi-
nation system is similar to that in humans. The structure of the female genome includes
the sex chromosome X. Individuals that have two female genomes are female (♀). The
structure of the male genome includes the sex chromosome Y. Individuals with one
female and one male genome are male (♂).
   Consider the simplest example of a L-E-HPS (Fig. 2), which includes a parent spe-
cies (P. lessonae) as well as hybrids which clonally transmit the female P. ridibundus
genome (XR). Hybrids are reproduced when crossed with the parental species individ-
uals: ♀XL(XR)×♂XLYL → ♀XL(XR) : ♂YL(XR); ♀XLXL×♂YL(XR) → ♀XL(XR).




Fig. 2. The occurrence of Pelophylax esculentus due to hybridization of Pelophylax lessonae
with Pelophylax ridibundus and reproduction of P. esculentus with the parental species in
L-E-HPS.

Different regions of the P. esculentus complex distribution are characterized by HPSs
of various compositions [2, 10, 11]. In some regions there are P. esculentus individuals
which simultaneously produce gametes (L) and (R). This phenomenon is called hybrid
amphispermy. Such individuals are referred as (L)(R).
   By crossing hybrids, which transmit genomes of the same parental species, there
may occur representatives of these species. Typically, these individuals die before sex-
ual maturity: ♀XL(XR)×♂YL(XR) → ♀XRXR → †.
   Different forms of frogs in a HPS differ in their vitality and fertility. The composition
of zygotes, tadpoles and frogs of the different ages in the same HPS may vary signifi-
cantly. Another feature of the hybridogenic complex of water frogs (not considered in
this study) is that in some regions there are not only diploid hybrids (i.e., having two
genomes), but hybrids with three (LLR or LRR) and even four (e.g. LLRR) genomes
[2, 15].
   Hemiclonal hybridization, the consequence of which is the occurrence of a HPS, is
observed not only for the water frogs, but also for some other species groups [1].
   HPSs of hybridogenic species complexes are a little-known category of biosystems.
Their study, part of which is a simulation of their transformations, should lead to im-
portant results. For example, the particular genome-selective elimination in interspe-
cific hybrids may open new opportunities in biotechnology and genetic medicine. This
allows, if necessary, the removal of unwanted fragments of the genome. Hemiclonal
inheritance supports all offspring from crossing of two individuals (each of which trans-
mits a clonal genome) to be genetically identical, i.e. clonal. The resulting organisms,
which inherit clonal genomes of two different parent species, may be useful for bio-
technology and agriculture. Study of the reason of the hemiclonal hybridization occur-
rence and the effects of this phenomenon on the evolution of hybridizing species are of
considerable theoretical interest. Although study of HPS transformations per se cannot
solve these problems, it allows one to better understand mechanisms of appearance and
maintenance of the stability of a biological system in which such amazing genetic phe-
nomena are possible.
   For further study, it is necessary to describe the variety of possible states of HPSs,
their dynamics, regularities and the conditions under which they are stable.
   Direct study of HPSs faces significant challenges. Determination of their composi-
tion and reproduction mechanisms is associated with a significant amount of fieldwork.
HPSs are relatively unusual and differ from well known biosystems such as biocenoses
and sexually or clonally reproducing populations. Processes of change in HPS are quite
extended temporally and may take decades. There is a high degree of variability in the
composition of different HPSs and some are unique objects whose occurrence is highly
unlikely. The study of only one unique HPS may not be enough for understanding the
regularities of their dynamics. Additionally, HPSs are complex systems containing
many stochastically interacting components.
   Some authors have used analytical modelling to describe the dynamics of a popula-
tion system of water frogs [14, etc.]. Analytical models are suitable for the study of
separate aspects of the HPS dynamics, however their usage faces significant challenges
due to fact that a HPS is governed by an interrelated set of stochastic processes. We
therefore consider simulation to be a more useful tool to study the general properties of
HPSs.
   Christiansen [4] simulated reproduction of P. esculentus. Christiansen’s model is
deterministic whereas our model allows modelling of random events in competition and
breeding of animals. The works of Bove et al. [3] and Quilodran et al. [12] are recent
studies based on the simulation of water frogs' HPSs. These studies consider the stabil-
ity of specific types of HPSs under certain conditions. Our work, in contrast, seeks to
analyze all possible stable states of a certain category of HPSs.
   Differences of our simulation from others are as follows. At the same time, we con-
sider typical for most populations demographic factors (non-competitive and competi-
tive mortality, differences in the reproduction probability, changes in viability with age)
and the unique features of HPS. Our work considers the entire space of possible states
of a certain category of HPS. Finally, we use population parameters, the estimates of
which were obtained during field studies of this category of HPS.

1.3    Scope of Modelling
Several authors of this work have previously built a deterministic discrete-time simu-
lation model [5]. It enabled investigators to determine the population composition for
a specified number of simulation steps, commencing with the current HPS composition.
The objective of the current work is identification of the set of stable states of a HPS
(consisting solely of diploid water frog individuals) using a stochastic simulation.
    The research was conducted in two stages. In the first stage a number of simulation
experiments were performed for various initial compositions of the model HPS. The
aim of this stage was to get insight into the possible end states for the HPSs. The second
stage analyzed and classified these end states. The dynamic typology of the system was
constructed based on the various initial compositions of the HPSs, the observed final
states, and the observed transitions from initial to final states.
    Dynamic typology is based not only on analysis of the observed object states, but,
above all, on a forecast of their future dynamics [9]. In this respect, dynamic typology
differs both from associative typology (identification of groups of objects, related to
one or more samples) and from analytical typology (partitioning of a set of objects into
groups depending on the state of their observable characteristics) [15].
    Hybrid frogs are reproduced differently in different regions [10, 14]. Before the
problem of the dynamical typology is solved in general, for all possible types of HPSs
and all known genetic forms of hybrids and their reproduction, it should be solved for
one selected region. Authors of the current paper have previously described the Siv-
erskyi Donets center of diversity of the Pelophylax esculentus complex located at the
eastern Ukraine [15]. It is characterized by a high diversity of HPSs and provides the
potential to test the adequacy of the model. The distribution of HPSs, including poly-
ploid hybrids (not discussed in this paper), is related to the flood plain of the river Siv-
erskyi Donets, and to its small tributaries and ponds located nearby within the Kharkiv
and Donetsk regions. Within the Mzha and the Uda river basins (right tributaries of the
Siverskyi Donets), HPSs consisting exclusively of diploids are widespread. The nature
of gamete production in diploids from the Siverskyi Donets and Mzha and the Uda river
basins is similar.
    We plan to describe the total variety of HPSs from the Siverskyi Donets center of
diversity of the Pelophylax esculentus complex. As the first step, we consider a variety
of HPSs located in the Mzha and the Uda river basins which consist of diploid repre-
sentatives.
2           Model Description and its Justification

2.1         The Cycle of the Model
A simulation model, as opposed to an analytical model, describes the process of the
state transformations during time, not the dependence of the future states of the system
on the current one. Therefore, for model development it is sufficient to describe the
algorithm of changes in the HPS sizes with time. In this model, time is divided into
discrete steps through which the model cycles. The cycle of the model corresponds to
the calendar year. It sets the sequence of transformations: αng a → βng a → γng a → δng a
→ ωng a . Here αng a , βng a , γng a , δng a and ωng a is a sequence of transformations of the
individual's groups of the genotype (g) and of the age (a) which correspond to the dif-
ferent stages of the annual cycle (Table 1).
   Of k genotypes, x is female and y is male one. At a certain stage, the females of
genotype f a and males of genotypes m a form a pair then descendants appear. The other
symbols are explained in the description of the input parameters of the model.

                               Table 1. A cycle of the model work.
                                                                          Changing the number of
    Symbol                Meaning                  Transformation
                                                                          individuals in the groups
                                                                                α g     ω     g
                     The initial numbers of        Start of cycle                t n a = t-1 n a
      α g          individuals in groups of a
        na
                 certain genotypes and age in Transition individuals                 β      g   α g
                                                                                         t n a = t n a-1
                           each cycle             to the next age
                The numbers of individuals in
      β g
        na        groups with allowance for
                their transition to the next age Non-competitive                γ         β g
                                                                                    tn a ≈ tn a × s a
                                                                                      g            g

                The numbers of individuals in          death
      γ g
        na       groups after non-competitive
                            mortality                                           δ      g   γ g        g
                                                   Immigration                      t n a = t n a+ t i a


      δ g       The numbers of individuals in                               Algorithm of the
        na
                groups after joining immigrants Competitive reduction calculation ω t ng a based on
                                                   in number of       δ g       g     g
                                                                       t n a , c a , d a and V (it is
                 The numbers of individuals in      individuals
                                                                           described below)
    ω g            groups after competitive
      na                                                                    Algorithm of the
                    mortality due to lack of
                           resources                                    calculation     P(ff a ,mm a' )
                                                  Creation of the
                                                                       based on ωnf a , lf a , ωnm a'
                                                   parental pairs
                                                                       and em a' (it is described
                  Number of pairs of a certain                             separately below)
P(ff a ,mm a' )
                          composition

                                                     Reproduction
                                                                            ω   g
                                                                             tn 0=        Σ(
                                                                                      P(ff a ,mm a' )×
                                                                             o (f ,m )×bf a ×wm a' )
                                                                              g f   m

    ω     g
        tn 0    The number of offspring (a=0)
                                                                     End of cycle
2.2    Input Parameters
Viability parameters. sga∈[0,1] — survival. This is proportion of individuals, which
happens to be saved in a result of non-competitive death. If а>=maxaga, then sga=0 where
max g
     a a stands for the maximum life span, and а is the age of individual.
    сg a ∈[0,1] — competitiveness factor. Let’s designate the probability of an individual
to survive over the competitive reduction as с´g a . Denote maxс´ as the maximum value
of this probability that is characteristic for representatives of the most competitive
groups (0<maxс´<=1). According to algorithm of the competitive reduction, the number
of the entire HPC along with the value of competitiveness affects a probability to sur-
vive for each individual. The value of the competitiveness factor
сg a =с´g a /maxс´ can be defined as one of the input parameters of the model.
    dg a , (dg a >=0) — demand. This is the number of resources required for individual of
a certain group. It’s magnitude is given for one cycle of the model.
    lf a ∈[0,1] — female loveliness. This is the female success rate in its search for a part-
ner for reproduction. That is set in the same way as the competitiveness factor is spec-
ified. Let’s denote the maximum value of the probability to find a partner as
max
     l´, 0<maxl´<=1). That is characteristic for the representatives of the most successful
groups of females. lf a =l´f a /maxl´. Denote the matag a as the age of sexual maturity. For
females with the age less than matag a (a<matag a ), lf a =0.
    Similarly, em a ∈[0,1] — male effectiveness. This is the male success rate in its search
for a partner for reproduction. Correspondingly, em a = e´m a / maxe´. For males, the age
of which (a<matag a ), em a =0.
    bf a ∈N — breed. This is the female fertility. That is specified by the number of eggs
produced by the females.
    wm a ∈[0,1] — wad. This is ability of the males to fertilize. It is measured by a fraction
of eggs that should be fertilized by a given male.
    og(ff,mm)∈[0,1] — offspring. This is the proportion of the descendants of the g-th
genotype of the female’s f-genotype with the male’s m-genotype. That is set for all
possible genotypes for entire pairs P(ff,mm).


Conditions for experiment. The following global settings should be set to run an ex-
periment or series of experiments.
   0 n a ∈N— number, which is the initial number of all groups.
        g

   t i a ∈N — inbound, which is number of individual-immigrants at a certain cycle (t)
      g

provided by the experiment script.
   V — volume — availability of resources that is amount of resources provided by
habitat. It is measured in the number of resource units.

2.3    The Algorithm of the Model
A cycle of the model work is described in the Table 1. The table describes variations in
the numbers of individual groups. However, in the model, calculations are made for
each individual and then summarized for the groups of individuals. For example, for a
non-competitive reduction of quantities γ t ng a ≈ β t ng a × sg a . This process is implemented
stochastically, and that explains usage of sign ≈, and not =. The pseudo-random number
generator defines the fate of each individual of the every genotype and every age. An
individual survives with probability sg a , and dies with probability 1-sg a .


Algorithm of the competitive mortality rate. Denote the δD = Σ(δtnga × dga),
ε g
 tn a = tn a × c a and D = Σ( tn a × d a ).
        δ g     g      ε            ε g    g
      δ            ω g δ g
   If D≤V, then t n a = t n a .
   If δD>V and εD=V, then ω t ng a ≈ ε t ng a .
   If δD>V and εD>V, then ω t ng a ≈ ε t ng a ×V/ εD.
   If δD>V and εD<V, then ω t ng a ≈ δ t ng a – (δ t ng a –ε t ng a )×(δD – V)/(V – εD).
   The fate of each individual of each genotype and each age is determined by the
pseudo-random number generator. They survive with probability ω t ng a / δ t ng a and die
with probability 1 – ω t ng a / δ t ng a .
   Such a reduction in the quantities satisfies the following conditions:
   1) Quantity before reduction is reduced up to the value that corresponds to resource
availability.
   2) Quantity of each group before the competitive reduction, δ t ng a , is reduced to a
quantity after competitive reduction, ω t ng a , in such a way that percent of individuals of
each group, which passed through the competitive reduction, is proportional to the com-
petitiveness of this group’s representative: ω t ng a / δ t ng a ~ cg a .


Algorithm of creation of the parental pairs. The quantity of the female’s ffa and the
male’s mma' pairs is determined in the following way:
   P(ff a ,mm a' ) ≈ ωnf a × (ωnf a ×lf a / Σ(ωnf a ×lf a )) × (ωnm a' ×em a' / Σ(ωnm a' )×em a' )).
   The total number of females from the group ff a is multiplied by the probability for
the females of the group to be selected by any male. This allows you to set the number
of the females in the group, who should find a partner in this iteration. The resulting
number is multiplied by the probability for the male of the group mm a' to be chosen by
any female. In this way one can identify the females from group ff a , which form a pair
with the male from group mm a' .
   The probability for the females from the group ff a to find a partner is determined by
the ratio of the product of quantity of this group ωnf a and the success factor of this group
of representatives in finding a partner lf a , to the sum of such products for all groups of
females.
   The probability for the males from the group mm a' to find a female is determined by
the ratio of the product of quantity of this group ωnm a' and the success factor of this
group of representatives in finding a partner em a' to the sum of such products for all
groups of males.
   The fate of each particular female and male is determined by the random number
generator in accordance with the probabilities, given the expected number of pairs of
P(ff a ,mm a' ).
3      Selection of Parameters

The adequacy of the model was tested on the processes observed in the all-diploid HPS
from the Siverskyi Donets center of diversity of the Pelophylax esculentus complex. In
total, eight different genomes can be transmitted in HPS of these water frogs (Table 2).
   In the Siverskyi Donets center of diversity of the Pelophylax esculentus complex
there are no mature P. lessonae, and therefore are no genomes of XL and YL. In this
paper, we have considered only diploid frogs. Therefore, in the simulation described
only eight genotypes were considered as shown in Table 3.

                      Table 2. The genomes considered in the model.

                        XL       YL       XR          YR
Genome                                                      (XL)        (YL)    (XR)   (YR)

Mode of inheritance              Recombinant                              Clonal

Species                 P. lessonae       P. ridibundus      P. lessonae       P. ridibundus

Sex                   Female    Male    Female     Male    Female       Male   Female Male


                      Table 3. The genotypes considered in this work.

                             Females                                Males

                               X X                                  X Y
P. ridibundus                  R R                                    R R

                       X
P. esculentus           R(XL), (XL)(XR)           X
                                                    R(YL), YR(XL), (XL)(YR), (YL)(XR).


The set of genomes is defined as G=[0,1]4×{0,1}2⊂R4. An individual is designated by
three parameters: age A∈N and two genomes (g1 and g2). Thus, a set of individuals is
A×G⊂R13.
   The set of HPS systems HPS=ФN+M⊂R13(N+M)consists of N + M frogs, where N is
the size of the spawning HPS, and M stands for the number of immature individuals.
   The following data input is used in the model:

— A description of the initial composition of the model HPS, and the scenario of mi-
  grants entering into it;
— A description of the parameters of the viabilities for all considered groups of indi-
  viduals that differ by age and by genomic composition;
— A description of the results of all possible crossings in the model (distribution of the
  probability of occurrence of offspring with a specific genomic composition from a
  crossing of different parents).
   A detailed justification of the chosen default settings is beyond the scope of this
work. It is done on the basis of study of the natural HPS, including above all the deter-
mination of the lifespan and the frog's growth rates using skeletochronology data, and
results from the estimation of the population size and composition of the natural HPS
by the mark-and-recapture method. Assumptions about the results of crossings (Table
4) were made based on the results of a study of the nature of gametogenesis of hybrid
frogs from the Siverskyi Donets center of diversity of the Pelophylax esculentus com-
plex [2]. The model assumes that resource consumption is proportional to the biomass
of the individual. Empty cells in the Table 4 correspond to the crossings of which the
offspring are non-viable. In cases when a cell, corresponding to a certain crossing, con-
tains 2, 3 or 4 genotype, the descendants that correspond to these genotypes appear with
equal probability (corresponding to ½, ⅓, or ¼).
   A preliminary research was carried out to determine the age structure of the model
HPS in which we determined the equilibrium age structure of the HPS P. ridibundus.
In further experiments, the age distribution of the initial composition of HPS corre-
sponds to that age structure. In those cases when we added the new specimens into HPS
at a certain step, they were considered as 4 years old.

                     Table 4. The proposed crossing results, assumed at modelling.

    Geno-
                 ♂♂             X Y
                                 R R       X
                                           R(YL)      Y
                                                      R(XL)           (XL)(YR)              (YL)(XR)
    types

                Gam-
        ♀♀                  X
                             R        Y
                                       R       (YL)       (XL)       (XL)    (YR)          (YL)    (XR)
                 etes

    X X          X       X X
        R R       R         R R XRYR XR(YL)           X
                                                      R(XL)      X
                                                                 R(XL)       X Y
                                                                              R R      X
                                                                                       R(YL)       X X
                                                                                                    R R

    X
    R(XL)        (XL)   X
                         R(XL) YR(XL)          —          —          —      (XL)(YR)       —      (XL)(XR)

                 (XL)   X
                         R(XL) YR(XL)          —          —          —      (XL)(YR)       —      (XL)(XR)
    X    X
( L)( R)
                 (XR)    X X
                            R R XRYR (YL)(XR) (XL)(XR) (XL)(XR) XRYR (YL)(XR) XRXR


4            Case Studies

4.1          The Space of Possible Outcomes of the HPS Transformations
To determine the set of possible end states of the water frogs' HPS, consisting of diploid
representatives, we chose a set of initial states, evenly distributed in the space of possi-
ble states. For 8-examined genotypes, the part of each (p♀XRXR, p♂XRYR etc.) varied
in increments of 1/8 from 0 to 7/8. Total 8! = 16,777,216 such combinations are possi-
ble. Of this amount, the combinations that satisfy the following conditions were se-
lected:
   p♀XRXR + p♂XRYR + p♀XR(XL) + p♂XR(YL) + p♂YR(XL) + p♀(XL)(XR) +
p♂(XL)(YR) + p♂(YL)(XR) = 1;
   p♀XRXR + p♀XR(XL) + p♀(XL)(XR) > 0;
   p♂XRYR + p♂XR(YL) + p♂YR(XL) + p♂(XL)(YR) + p♂(YL)(XR) > 0.
   Naturally, it does not make sense to consider the combinations in which the total
proportion of genotypes involving into the HPS is not equal to unity, as well as those
in which there are no male or female. The total number of genotypes that satisfy the
above conditions equals 5895. Ten simulations of 500 steps were conducted for the
5895 starting points. The collection of all observed outcomes of the simulation was
divided into types, depending on what kind of genotypes presented in the model HPS
over 500 steps. For 4778 initial states, all 10 iterations led to any one outcome. In 1117
cases, outcomes were variable. For these 1117 initial states, 10 more runs were carried
out; thus, the total number of simulations was 70,120.
   To determine the states, which the model HPS can move to, we examined the inter-
vals between step 100 and 200, step 200 and 300, 300 and 400 as well as 400 and 500.
The total number of intervals were 162,580; its number is less than possible one due to
the model HPS, which dies at any stage of the simulation.
   To divide the obtained set of 70,120 final states of individual simulations into
groups, we used the analysis of this aggregation by the method of principal components.
The first and the second major components make it possible to divide the HPS into 6
groups (Fig. 3. A):
— Extinction — 37,946 runs ended in HPS extinction;
— R-E-HPS — 22,204 runs ended in different HPSs R-E-types;
— E-HPS-type I — 5884 results; such HPS include ♀(XL)(XR) and ♂(XL)(YR);
— E-HPS-type II — 3892 results corresponding to another possible type of the E-HPS,
   which include ♀(XL)(XR) and ♂(YL)(XR);
— Extincting — 124 results located nearby of the point Extinction. They continue to
   end up in extinction. The structure of these finals can be further divided into 13
   types;
— R-population — 70 results, consisting only of the population P. ridibundus.
   The R-E-HPS group, which corresponds to 22,204 results by its composition, can be
divided into three parts. Their relative positions can be seen in the plane of the first and
the third main component (Fig. 3. B):
— Stable R-E-HPS-type I — 18,852;
— Stable R-E-HPS-type II — 2,847 results;
— Indifferent R-E-HPS — 505 results in which P. ridibundus is present as well as
P. esculentus, which transmits both genome P. lessonae and female genome P. ridi-
bundus.
Fig. 3. Ordination of the results of 70,120 simulations. A. Ordination of the results on
the plane of the first two principal components. B. Ordination of the results on the plane
of the first and third principal components.

To determine the types of stability for the dynamic types of the allocated HPS, the 162
580 pairs of initial state and its outcome were analyzed. The observed stated were clas-
sified by the types of stability as shown in Fig. 4.
    The observed states are divided into three groups (Equilibrium states, Transient
states and Attractive states) as shown in Fig.6, though this division is rather subjective.
The Indifferent equilibrium (II) can be considered as a transient state as well. Two ver-
sions of the transient states are associated with the directed transitions to other states.
Extincting state (IV) is associated with the transition to the extinction state (VI), and
transforming stage (III) with the transition to other states.




   Fig. 4. Types of biosystems stability, observed in experiments with the simulation model.

Six basins of sustainability exist in the space of possible states of HPS that consists of
diploid individuals, and this corresponds to specific features of the Siverskyi Donets
center of diversity of the Pelophylax esculentus complex by their character of the ge-
nome transmission. One of them complies with the population of parental species. Two
are E-HPS, i.e. spawning population of HPS consisting exclusively of P. esculentus.
When crossing P. esculentus with the hybrid amphispermy, the offspring of parental
species appearing in such systems dies before the age of maturity.
   E-HPS-type I:
   ♀(XL)(XR)×♂(XL)(YR)→XLXL:(XL)(XR):(XL)(YR):XRYR→♀(XL)(XR) : ♂(XL)(YR);
   E-HPS-type II:
   ♀(XL)(XR)×♂(YL)(XR)→XLYL:(XL)(XR):(YL)(XR):XRXR→♀(XL)(XR) : ♂(YL)(XR).
   There are two more basins of stability corresponding to the R-E-HPS containing
P. ridibundus and P. esculentus, and the last basin is associated with extinction of the
HPS.


5       Outcome Interpretation

When assessing the results, one must be aware that the result of the simulation is not
proof of a hypothesis. At the same time, using simulation models as an exploratory tool
has a significant advantage over unformulated conceptual models. We do not have the
sufficient empirical data to describe precisely the processes occurring in the natural
HPS. The lack of empirical data is offset by a set of presumptions [13] and hypotheses.
Simulation allows us to derive consequences arising from the set of initial assumptions.
   These sets of consequences may or may not contradict the observable empirical pic-
ture. An observed contradiction is the basis for rejecting a set of initial assumptions or
adjustments. Agreement is not an evidence for the initial assumptions, but can be seen
as an argument in their favor, in other words as corroboration, not as proof (Fig. 5).




Fig. 5. Using a simulation model to test hypotheses about the mechanisms of HPS functioning.
The distribution of the outcomes of the model HPS agrees with the empirical data, and it is ob-
tained in the case when variant B of the variable part of the model is selected. It does not prove
variant B to be true, but it allows to discard variant A in its favor.
Fig 5 shows how the simulation results can be used to select between different versions
of the assumptions in the absence of empirical data. Model predictions regarding the
expected diversity of HPS conditions, depends on assumptions, on which the variant
part of the model was constructed. The predicted (modeled) expected diversity of HPS
states can be compared to empirically observed diversity of HPS states. The results of
this comparison can be seen as the arguments in favor of such initial assumptions that
yielded close to the empirical distribution of states HPS.
    Simulation results, among other things, stimulate the collection of the empirical data.
Comparison of predictions, obtained by modelling, and the available fragmented data
about the composition of the natural HPS and character of P. esculentus gametogenesis
shows that currently there is an absence of E-type HPS in the studied region, although
that there are grounds for thinking that such a system probably existed in 1995 [6].
Studying the composition of tadpoles in this pond shows that there existed E-HPS-type
II [8]. In fact, the majority of observed systems are in a transient state.
    The authors believe that this research is only a starting step in design of the dynamic
typology of the water frogs’ HPS. The usefulness of the simulation is not limited to fact
that two stable states, previously not known to the authors, can theoretically be found:
R-E-HPS-type II and E-HPS-type II. Results of simulations allow developing a pro-
gram for further research. The simulation results define the data collection program for
testing the adequacy of the results. The currently available data do not contradict the
model findings. Future work will include the study of unstable states of HPS that are
observed in the natural environment [16], as well as the extension of the model by in-
corporating hybrids with three genomes, which are specific for the Siverskyi Donets
center of diversity of the Pelophylax esculentus complex.


6      Discussion and Conclusion

We study the unusual category of biosystems, Hemiclonal Population Systems (HPS),
through the example of the hybridogenic complex of water frogs, the Pelophylax escu-
lentus complex. The unusual method of reproduction of interspecific hybrids within
HPS results in their unusual features, which need more investigation. An important
method for studying such systems is computer simulation.
   A simulation model of the water frogs’ HPS has been presented. The model inputs
are the parameters describing the comparative vitality of various genetic forms of frogs,
the results of their probable crossings, as well as the experimental conditions such as
the capacity of the environment, the initial composition of the model HPS, and a sce-
nario for the introduction of migrated frogs into the model HPS. Evaluation of the com-
parative vitality and the crossing results are defined in accordance with the results of
population-environmental research in the Siverskyi Donets center of diversity of the
Pelophylax esculentus complex. Since direct empirical data was lacking, the values of
the parameters used have been estimations that authors put forward on the basis of a
study of this region.
   Repeated runs of the model yield a probability distribution of outcomes of various
HPS transformations according to their initial states and experimenter-determined pa-
rameters. We analyzed the processes in all-diploid HPS from the Siverskyi Donets cen-
ter of diversity of the Pelophylax esculentus complex, which is characteristic for the
Mzha and the Uda river basins. When sampling evenly among the starting point of runs
in the space of possible states of such systems, a set of outcomes has been established.
Only six stable states have been found for such systems, one of which corresponds to
the population of the parental species, two coincide with the HPS consisting of repre-
sentatives of the parental species and hybrids, two correspond to the HPS consisting
solely of hybrid individuals, and one state corresponds to the extinction of frogs in the
simulated habitat.


References
 1. Abbott, J.K., Morrow, E.H. Obtaining snapshots of genetic variation using hemiclonal anal-
    ysis. Trends in Ecology & Evolution 26(7), 359–368 (2011).
 2. Biriuk, O.V., Shabanov, D.A., Korshunov, A.V., Borkin, L.J., Lada, G.A., Pasynkova, R.A.,
    Rosanov, J.M., Litvinchuk, S.N: Gamete production patterns and mating systems in water
    frogs of the hybridogenetic Pelophylax esculentus complex in northeastern Ukraine. J. of
    Zool. Syst. & Evol. Res., 54(2), doi:10.1111/jzs.12132 (2016).
 3. Bove, P., Milazzo, P., Barbuti, R.: The role of deleterious mutations in the stability of hy-
    bridogenetic water frog complexes. BMC Evolutionary Biology, 14:107, doi:10.1186/1471-
    2148-14-107 (2014).
 4. Christiansen, D.G. Gamete types, sex determination and stable equilibria of all-hybrid pop-
    ulations of diploid and triploid edible frogs (Pelophylax esculentus). BMC Evolutionary Bi-
    ology, 9, doi:10.1186/1471-2148-9-135 (2009).
 5. Kravchenko, M.O., Shabanov, D.A., Vladimirova, M.V., Zholtkevych, G.M.: Investigation
    of the stability of hemiclonal population systems of water frogs hybridogenetic complex by
    the means of simulation modeling. The Journal of Dnipropetrovsk National University. Bi-
    ology. Ecology (in Russian) (2011).
 6. Lada, G.A.: On the necessity of preserving the unique "pure" diploid populations of the
    edible frog (Rana esculenta Linnaeus, 1758) in the Belgorod and the Kharkiv region. In:
    Problems of protection and rational use of natural ecosystems and biological resources.
    Penza (in Russian) (1998).
 7. Lenhard, J.: Computer Simulation: The cooperation between experimenting and modelling.
    Philosophy of Science, 74(2), 176-194 (2007).
 8. Makaryan, R.M., Biryuk, O.V., Korshunov, O.V., Kravchenko, M.O., Meleshko O.V.,
    Trokhimchuk R.R., Shabanov D.A.: Composition of water frogs (Pelophylax esculentus
    complex) tadpoles in the Is’kov Pond (NPP "Gomilshanskie Forests"). In: Status and biodi-
    versity of the ecosystems of the Shatsk NNP & other reserved areas. Lviv (in Ukrainian)
    (2016).
 9. Melechov, I.S.: Dynamic typology of forests. Forestry (in Russian) (1968).
10. Mezhzherin, S.V., Morozov-Leonov, S.Yu., Rostovskaya, O.V., Shabanov, D.A., Sobo-
    lenko, L.Yu.: The ploidy and genetic structure of hybrid population of water frogs Pelophy-
    lax esculentus complex (Amphibia, Ranidae) of Ukraine fauna. Cytology and Genetics,
    44(4), 212–216 (2010).
11. Plötner, J.: Die westpaläarktichen Wasserfrösche. Bielefeld: Laurenti Verlag (2005) .
12. Quilodran, C.S., Montoya-Burgos, J.I., Currat, M.: Modelling interspecific hybridization
    with genome exclusion to identify conservation actions: the case of native and invasive Pel-
    ophylax waterfrogs. Evolutionary Applications, 8 199–210 (2015).
13. Rasnitzyn, A.P.: Theoretical Foundations of Evolutionary Biology. In: Introduction to pal-
    aeoentomology. Moscow: KMK (in Russian) (2008).
14. Reyer, H.-U., Wälti, M.-O., Bättig, I., Altwegg, R., Hellriegel, B.: Low proportions of re-
    producing hemiclonal females increase the stability of a sexual parasite–host system (Rana
    esculenta, R. lessonae). Journal of Animal Ecology, 73, 1089–1101 (2004).
15. Shabanov, D.A.: Evolutionary ecology of population hybridogenic complex of water frogs
    (Pelophylax esculentus complex) Left-Bank Ukraine steppe: Thesis for the Degree of Doctor
    of biological sciences, spec. 03.00.16 ecology. Dnipropetrovsk (in Ukrainian) (2015).
16. Shabanov, D., Usova, O., Kravchenko, M., Biriuk, O., Leonov, A., Korshunov, O., Mair,
    Q., Meleshko, O., Newman, J., Vladymyrova, M., Zholtkevych, G.: Sustainable coexistence
    of the parental species and hemiclonal interspecific hybrids is provided by the variety of
    ontogenetic strategies. Herpetological Facts Journal, 2, 35–43 (2015).

</pre>