Alejandro Lopez Rincon

Carmina A Perez Romero
email: perezr.carmina@gmail.com

Lucero Mendoza Maldonado

Juan I Menchaca

Eric Claassen
email: prof.eric.claassen@gmail.com

Johan Garssen
email: johan.garssen@danone.com

Aletta D Kraneveld
email: a.d.kraneveld@uu.nl

Alberto Tonda
email: alberto.tonda@inrae.fr

Design of Specific Primer Sets for SARS-CoV-2 Variants Using Evolutionary Algorithms

Keywords: Applied computing → Genomics; Bioinformatics primers, genetic algorithm, SARS-CoV-2, B.1.1.7

Primer sets are short DNA sequences of 18-22 base pairs, that can be used to verify the presence of a virus, and designed to attach to a specific part of a viral DNA. Designing a primer set requires choosing a region of DNA, avoiding the possibility of hybridization to a similar sequence, as well as considering its GC content and 𝑇 𝑚 (melting temperature). Coronaviruses, such as SARS-CoV-2, have a considerably large genome (around 30 thousand nucleotides) when compared to other viruses. With the rapid rise and spread of SARS-CoV-2 variants, it has become a priority to breach our lack of specific primers available for diagnosis of this new variants. Here, we propose an evolutionary-based approach to primer design, able to rapidly deliver a high-quality primer set for a target sequence of the virus variant. Starting from viral sequences collected from open repositories, the proposed approach is proven able to uncover a specific primer set for the B.1.1.7 SARS-CoV-2 variant. Only recently identified, B.1.1.7 is already considered potentially dangerous, as it presents a considerably higher transmissibility when compared to other variants.

INTRODUCTION

As the pandemic of SARS-CoV-2 continues to affects the globe, researchers and public health teams around the world monitor the virus for acquired mutations that may lead to higher risks of developing COVID-19 or vaccine resistance. Prompt and widespread diffusion of information related to viral threats plays a critical role in research and mitigation of outbreaks [START_REF] Mark B Van Der Waal | Blockchain-facilitated sharing to advance outbreak R&D[END_REF]. This is especially true when viruses mutate rapidly under clinical, therapeutic or vaccine pressure.

Although SARS-CoV-2 acquires 1-2 mutations every month [START_REF] Duchene | Temporal signal and the phylodynamic threshold of SARS-CoV-2[END_REF][START_REF] Lu | Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding[END_REF], a new variant has been recently reported in the UK as a Variant of Concern (VOC) 202012/01 [START_REF] Chand | Investigation of novel SARS-COV-2 variant: Variant of Concern 202012/01[END_REF][START_REF]Rapid increase of a SARS-CoV-2 variant with multiple spike protein mutations observed in the United Kingdom[END_REF], from the B.1.1.7 lineage [START_REF] Arambaut | Preliminary genomic characterisation of an emergent SARS-CoV-2 lineage in the UK defined by a novel set of spike mutations[END_REF][START_REF] Chand | Investigation of novel SARS-COV-2 variant: Variant of Concern 202012/01[END_REF], Nextstrain clade 20B [START_REF] Hadfield | Nextstrain: real-time tracking of pathogen evolution[END_REF], or clade GR from GISAID (Global Initiative on Sharing All Influenza Data) [START_REF] Shu | GISAID: Global initiative on sharing all influenza data-from vision to reality[END_REF]. The B.1.1.7 variant presents 14 non-synonymous mutations, 6 synonymous mutations, and 3 deletions. The multiple mutations present in the viral RNA encoding for the spike protein (S) are of most concern, such as; deletion △69-70, deletion △144, N501Y, A570D, D614G, P681H, T716I, S982A, D1118H [START_REF] Arambaut | Preliminary genomic characterisation of an emergent SARS-CoV-2 lineage in the UK defined by a novel set of spike mutations[END_REF][START_REF]Rapid increase of a SARS-CoV-2 variant with multiple spike protein mutations observed in the United Kingdom[END_REF]. The SARS-CoV-2 S protein mutation N501Y alters the protein interactions involved in receptor binding domain. The N501Y mutation has been shown to enhance affinity with the host cells ACE2 receptor [START_REF] Arambaut | Preliminary genomic characterisation of an emergent SARS-CoV-2 lineage in the UK defined by a novel set of spike mutations[END_REF][START_REF] Tyler N Starr | Deep mutational scanning of SARS-CoV-2 receptor binding domain reveals constraints on folding and ACE2 binding[END_REF] and to be more infectious in mice [START_REF] Gu | Adaptation of SARS-CoV-2 in BALB/c mice for testing vaccine efficacy[END_REF].

The mutation P681H is situated in key residues of the S protein needed for viral transmission and viral entry in lung cells [START_REF] Arambaut | Preliminary genomic characterisation of an emergent SARS-CoV-2 lineage in the UK defined by a novel set of spike mutations[END_REF][START_REF] Hoffmann | A multibasic cleavage site in the spike protein of SARS-CoV-2 is essential for infection of human lung cells[END_REF].

The rate of transmission of the B.1.1.7 SARS-CoV-2 variant has been estimated to be 56%-70% higher, and its reproductive number (Rt) seems to be up to 0.4 higher [START_REF] Chand | Investigation of novel SARS-COV-2 variant: Variant of Concern 202012/01[END_REF][START_REF] Davies | Estimated transmissibility and severity of novel SARS-CoV-2 Variant of Concern 202012/01 in England[END_REF]. Since it first identification in southern England, the presence of the B.1.1.7 variant has been rapidly increasing in the UK, displacing other SARS-CoV-2 variants [START_REF] Arambaut | Preliminary genomic characterisation of an emergent SARS-CoV-2 lineage in the UK defined by a novel set of spike mutations[END_REF]. Recently, Public Health England reported that infection with the B.1.1.7 variant can lead to an increase in risk of death of 1.65 when compared to infection with other non-VOC viruses [START_REF] Horby | NERVTAG note on B.1.1.7 severity[END_REF]. Furthermore, 11 samples of the B.1.1.7 variant have been found to carry the E484K mutation in the S protein, which has been shown to moderately enhance binding affinity of the ACE2 receptor [START_REF] Tyler N Starr | Deep mutational scanning of SARS-CoV-2 receptor binding domain reveals constraints on folding and ACE2 binding[END_REF] and to significantly reduce convalescent serum neutralization [START_REF] Allison J Greaney | Comprehensive mapping of mutations to the SARS-CoV-2 receptor-binding domain that affect recognition by polyclonal human serum antibodies[END_REF][START_REF] Allison J Greaney | Complete mapping of mutations to the SARS-CoV-2 spike receptor-binding domain that escape antibody recognition[END_REF], raising concerns about the current monoclonal/polyclonal antibodies used for treating the COVID-19 disease, as well as the efficacy of new mRNA vaccines in combating the spread of this new SARS-CoV-2 variants [START_REF] Gard Nelson | Molecular dynamic simulation reveals E484K mutation enhances spike RBD-ACE2 affinity and the combination of E484K, K417N and N501Y mutations (501Y.V2 variant) induces conformational change greater than N501Y mutant alone, potentially resulting in an esca[END_REF][START_REF] Wang | mRNA vaccineelicited antibodies to SARS-CoV-2 and circulating variants[END_REF]. B.1.1.7 and other N501Y-carrying SARS-CoV-2 variants have also been identified in different parts of Europe, South Africa, Brazil, Japan, Australia, USA, and Egypt [START_REF] Alm | Geographical and temporal distribution of SARS-CoV-2 clades in the WHO European Region, January to June 2020[END_REF][START_REF] Elghazaly | Laboratory based Retrospective Study to determine the start of SARS-CoV-2 in Patients with Severe Acute Respiratory Illness in Egypt at El-Demerdash tertiary hospitals[END_REF][START_REF] Hadfield | Nextstrain: real-time tracking of pathogen evolution[END_REF][START_REF] Shu | GISAID: Global initiative on sharing all influenza data-from vision to reality[END_REF].

Several molecular kits have been proposed and developed to diagnose SARS-CoV-2 infections. Most kits rely on the amplification of one or several genes of SARS-Cov-2 by real-time reverse transcriptase-polymerase chain reaction (qRT-PCR) [START_REF] Afzal | Molecular diagnostic technologies for COVID-19: Limitations and challenges[END_REF][START_REF]Molecular assays to diagnose COVID-19: summary table of available protocols[END_REF], using primers. Primers are short RNA sequences, used to detect the presence of a specific virus in a host organism. Good qRT-PCR primer candidates need to possess several qualities, including having a PCR product size or amplicon between 100 -200 bps, with a melting temperature (𝑇 𝑚 ) close to 60°C, a presence of C and G bases representing 40-60% of the sequence, and of a typical length of 18-22 nucleotides.

Public Health England was able to identify the surge of the B.1.1.7 SARS-CoV-2 variant only through the increase in S-gene target failure (negative results) from the otherwise positive target genes (N, ORF1ab) in their three target gene assay [START_REF] Chand | Investigation of novel SARS-COV-2 variant: Variant of Concern 202012/01[END_REF]. In other words, the primers designed for SARS-CoV-2 are not suited to identify its B.1.1.7 variant, due to fundamental differences in the RNA structure of the two virus strains. In order to detect and contain the potentially dangerous B.1.1.7, new primers need to be developed.

Primer design is traditionally a long process, where domain experts create a canonical sequence for the target virus, analyze the virus' structure, identify promising areas that might be unique to the strain, and create candidates that are then tested first in-silico, and then in the lab. The emerging adoption of machine learning techniques in the health sector, together with the free access to an unprecedented amount of data, has recently led to the first appearance of semi-and fully automated solutions for primer design, based on deep learning techniques [START_REF] Lopez-Rincon | Design of Specific Primer Set for Detection of B.1.1.7 SARS-CoV-2 Variant using Deep Learning[END_REF][START_REF] Lopez-Rincon | Classification and specific primer design for accurate detection of SARS-CoV-2 using deep learning[END_REF]. While effective, however, these approaches present several limitations, such as the inability of handling constraints on the candidate primers, and the training time required.

In this paper, we propose a novel, evolutionary-based approach to primer generation. An evolutionary algorithm (EA) [START_REF] De | Evolutionary Computation: A Unified Approach[END_REF] is applied to the task of finding meaningful sub-sequences of RNA that could represent good primer candidates. The fitness function is able to include constraints ranging from achieving a melting temperature in a given range, to presenting a sufficient quantity of specific bases.

Tested on the B.1.1.7 SARS-CoV-2 variant, the approach proves to be much faster than previous machine learning techniques, delivering results of the same quality.

The rest of the paper is organized as follows: Section 2 provides minimal background information on primers and their design. Section 3 outlines the proposed approach to primer design, based on an evolutionary algorithm. Section 4 summarizes the experimental evaluation, while Section 5 discusses the results and Section 6 concludes the paper.

BACKGROUND

This section briefly introduces the notion of primers, and provides an overview of the techniques used to design them.

Primers

A primer or oligonucleotide is a short nucleic acid sequence that helps to start the DNA synthesis in a Polymerase Chain Reaction (PCR) assay. The presence of primers is necessary because the DNA polymerase can only attach new nucleotides to an existing strand of nucleotides to make copies of a DNA fragment 1 . Primers typically have a length of 18-22 nucleotides and are essential components of the procedure, as they determine the specificity of PCR [START_REF] Van Pelt-Verkuil | PCR Primers[END_REF]. In the case of SARS-CoV-2, the viral RNA will be first transformed into single-strand complementary DNA (cDNA), and then the section of interest we will amplified. The forward and reverse primer attach in different directions of the targeted section in the cDNA sequence. If the cDNA sequence does not contain the forward and reverse primer sequences, then they will not attach and consequently won't be amplified, see Fig. 1 for a summary of the procedure.

Primer design is a crucial step of the process, because of its relationship with the success and quality of PCR analyses. Different factors can affect primers' efficiency, including; dimer formation (high self-complementary formation), stem loop interference, GC content (presence of C and G bases in the DNA sequence), high 3' stability (presence of C/G pairs at the 3' end) and extreme melting temperatures. While different computer programs have been proposed to facilitate the design process, it is still necessary to develop new algorithms to select the best PCR primers and avoid errors when using them in the laboratory [START_REF] Li | Primer design for RT-PCR[END_REF][START_REF] Rodríguez | Design of primers and probes for quantitative real-time PCR methods[END_REF].

Approaches to primer generation

The traditional approach to primer generation is to first use a canonical sample of the target variant, a summary of all samples available, and then analyze the differences with respect to the original virus. The differences, mainly mutations, are going to be possible locations of interest, and candidate primers can be designed by hand around the mutation. In the case of SARS-CoV-2, all the mutations are referenced to the canonical sequence 𝑁𝐶_045512.2 [START_REF] Wu | A new coronavirus associated with human respiratory disease in China[END_REF]. In addition, several parameters need to be verified such as 𝑇 𝑚 close to 60°, %GC content close to 50%, self complementarity, and high 3' stability. Then, the designed primer has to be checked as a unique sub-sequence, not appearing more than once inside the same sample. This is traditionally carried out with the BLAST software, from the National Center for Biotechnology Information (NCBI) [START_REF] Grzegorz M Boratyn | BLAST: a more efficient report with usability improvements[END_REF]. Nevertheless, if a sequence is not available in BLAST it cannot be identified using this method, as is the case for many of the new samples recently sequenced from the SARS-CoV-2 virus and its variants.

More recently, deep learning (DL) techniques for primer design have been applied to SARS-CoV-2 [START_REF] Lopez-Rincon | Classification and specific primer design for accurate detection of SARS-CoV-2 using deep learning[END_REF]. The DL approach is based on a Convolutional Neural Network (CNN) that is trained to separate the target virus' samples from the rest. Then, the filters of the CNN are analyzed to extract sequences of 21 base pairs, that are later tested for the suitability as primer candidates using the software Primer3Plus [START_REF] Untergasser | Primer3Plus, an enhanced web interface to Primer3[END_REF]. The same technique has been also applied to finding primers for the B.1.1.7 (British), the B.1.351 (South African), and the P.1 (Brazilian) SARS-CoV-2 variants. Around 16 hours are needed to find a suitable primer for each variant. Furthermore, most of the sequences generated by the DL approach are discarded due to undesirable properties (low presence of GC bases, high/low melting temperature, etc.) [START_REF] Lopez-Rincon | Design of Specific Primer Set for Detection of B.1.1.7 SARS-CoV-2 Variant using Deep Learning[END_REF].

PROPOSED APPROACH

In this work, we present a novel evolutionary-based approach to primer generation. As each primer candidate can be seen as a subsequence of a SARS-CoV-2 RNA sequence, our aim is to efficiently explore the search space of all possible sub-sequences among all available samples of a target virus strain. The ideal primer candidate would be present in all samples of the target strain, and not present in any other virus, be it a variant of the target, or another virus entirely. Furthermore, the candidate should possess a high number of C and G bases in its genome, a number of N values (indicating sequencing errors) as small as possible, and a temperature of melting 𝑇 𝑚 as close as possible to 60°.

Individual representation

An individual in our approach is a primer candidate, thus a subsequence of a given virus sample. Assuming we have 𝑁 samples of the target virus, and each sample is a sequence of 𝐿 bps, an individual will be described by two integers, sample number 𝑘 in the training dataset, and position number 𝑝 inside the sample. The candidate primer will then be represented by the 21 bps in positions {𝑝, 𝑝 + 1, ..., 𝑝 + 20} of sample 𝑠 𝑘 , see Figure 2. While primers can be of any length between 18 and 22 bps, we selected 21 to be able to compare our final results with the primers produced by DL techniques [START_REF] Lopez-Rincon | Design of Specific Primer Set for Detection of B.1.1.7 SARS-CoV-2 Variant using Deep Learning[END_REF][START_REF] Lopez-Rincon | Classification and specific primer design for accurate detection of SARS-CoV-2 using deep learning[END_REF], that are of length 21.

Population initialization

The population is initialized with random individuals. More specifically, the 𝑖-th individual will be characterized by two integers, 𝑘 𝑖 drawn with uniform probability in (1, 𝑁 ) and 𝑝 𝑖 drawn with uniform probability from (1, 𝐿 -21).

Fitness Function

The fitness function is a weighted sum, attempting to take into account all the criteria for a good primer candidate. The first term of the sum, 𝐹 , is evaluating the presence of the sequence selected as candidate primer 𝐼 inside training samples labeled with the target virus strain, and its absence from samples with a different label. In formal terms:

P (𝐼 ) = 𝑇 𝑖=0 𝑃 (𝐼, 𝑠 𝑖 ) ( 1
)
where 𝑇 is the number of samples in the training set, 𝑠 𝑖 is the 𝑖-th sample in the training set, and function 𝑃 is defined as:

𝑃 (𝐼, 𝑠 𝑖 ) = 0, if 𝐼 is found inside 𝑠 𝑖 and 𝐿(𝑠 𝑖 ) == 𝐿(𝑠 𝑘 ) 1, otherwise. ( 2
)
where 𝐿(𝑠) returns the class label of sample 𝑠. In other words, 𝑃 (𝐼, 𝑠 𝑖 ) equals 1 if sequence 𝐼 is found inside a sample with the same class label as sample 𝑠 𝑘 , the origin of sequence 𝐼 . So, if the candidate primer 𝐼 is found inside a sample that does not belong to the target class, or is not found in a sample that belongs to the target class, the solution is penalized.

The second term of the weighted sum takes into account the GC content of the candidate primer, or in other words, the presence of bases G and C:

C(𝐼 ) = 0.5 - 21 𝑖=0 𝐶 (𝐼 (𝑖)) 21 where 𝐶 (𝑏) = 1, if base 𝑏 is C or G 0, otherwise. (3)
where 𝐼 (𝑖) represents the base in position 𝑖 inside sequence 𝐼 . The following element of the weighted sum is N , defined as: where C(𝐼 ) is the content of C and G bases in sequence 𝐼 , as described in Equation 3, [𝑁 𝑎+] is the molar sodium concentration, and 𝑙 (𝐼 ) is the length of sequence 𝐼 , in bps. We used the value of [𝑁 𝑎+] = 0.2 as described in [START_REF] Untergasser | Primer3Plus, an enhanced web interface to Primer3[END_REF], while 𝑙 (𝐼 ) = 21 by design. The term taking into account 𝑇 𝑚 will then be:

N (I) =
T (𝐼 ) = |60 -𝑇 𝑚 (𝐼 )| (6)
To summarize, the final weighted sum will be:

F (𝐼 ) = 𝑤 𝑝 • P (𝐼 ) + 𝑤 𝑐 • C(𝐼 ) + 𝑤 𝑛 • N (𝐼 ) + 𝑤 𝑡 • T (𝐼 ) (7)
with 𝑤 𝑝 , 𝑤 𝑐 , 𝑤 𝑛 , 𝑤 𝑡 representing the weights associated to each term. Function F is to be minimized.

Selection

Individual selection for reproduction is performed through a tournament selection of size 𝜏, where 𝜏 individuals are randomly drawn from the population with uniform probability, their fitness is compared, and the individual with the best fitness is ultimately chosen.

Evolutionary operators

Mutations operators can act on the two integer values 𝑘, representing the index of a sample in the training set, and 𝑝, representing a position inside sample 𝑘's genome. Mutations acting on 𝑝 draw from a multinomial distribution where the probabilities are calculated from a normal distribution, following the idea that small displacements of the beginning of a primer left or right in the genome might provide small changes in the fitness function. Mutations targeting 𝑘, on the other hand, draw uniformly from all available sample indexes: most genomes of a virus, resulting from a sequencing process, will be almost aligned with each other, with small differences resulting from additions or deletions in the cDNA; for this reason, we believe that the principle of locality will be preserved when creating a primer candidate from position 𝑝 of two different samples 𝑘 and 𝑘 ′ .

The proposed crossover is a one-point crossover

Replacement

Individual replacement is carried out through a (𝜇, 𝜆) scheme, where the entire population is replaced at each generation.

EXPERIMENTAL EVALUATION

The proposed approach is validated on real-world data, using samples collected from the GISAID repository. 10,712 SARS-CoV-2 sequences are downloaded on December 23rd, 2020. After removing repeated sequences, we obtain a total of 2,104 sequences labeled as B.1.1.7, and 6,819 sequences from other variants, for a total of 8,923 samples. B.1.1.7 variants are assigned class label 1, while all the remaining samples are labeled as class 0. Then, we divided the data into 8,030 samples for training and 893 for testing. The algorithm described in Section 3 has been implemented in C#, and is available, along with all the data, in a GitHub open repository2 . After a few preliminary trials to tune the EA's parameters, population size is set to 𝜇 = 100, offspring size to 𝜆 = 100, and a termination condition based on a maximum number of generations, 20, probability of mutation to 𝑝 𝑚 = 0.15 and probability of crossover to 𝑝 𝑐 = 0.85, tournament selection size 𝜏 = 2, and a (𝜇, 𝜆) scheme, where the whole population is replaced at every generation. For the fitness function, the values for the weights are set as 𝑤 𝑝 = 1.0, 𝑤 𝑐 = 100.0, 𝑤 𝑛 = 1000.0, 𝑤 𝑡 = 1.0, to provide the same numerical importance to different parts of the evaluation, with the exception of N (𝐼 ), as primers containing N symbols are unlikely to be acceptable. In principle, the fitness function could also be structured as a lexicographic comparison, first evaluating the presence of Ns in the sequence, and then proceeding with the rest of the evaluation. Using a combination of lexicographic acceptability and multi-objective evaluations could also be beneficial. These possibilities will be analyzed in future works.

In order to assess the stability of the results, the proposed approach is run 20 times, with each single run lasting around 62 minutes with 5 threads on a 64-bit Windows 10 laptop with Intel Xeon E-2186M. Thread parallelization is exploited to evaluate several individuals at the same time. The best individual of each run is then stored, to perform a later comparison with primers designed through other methods, discussed in Section 5.

The 20 runs of the proposed approach return different candidate primers, as shown in Table 1. Each of the resulting solutions is then simulated in the canonical version of B.1.1.7 in the GISAID sample 𝐸𝑃𝐼 _𝐼𝑆𝐿_601443 using Primer3Plus [START_REF] Untergasser | Primer3Plus, an enhanced web interface to Primer3[END_REF], to compute the reverse primer and perform a final assessment of its suitability as a primer. This last process is more computationally expensive than a simple evaluation of a fitness function, and it can last up to 20 minutes (Primer3Plus) depending of the software used to crosscheck.

From the 20 best individuals obtained, only 6 solutions can be used as forward primers. Further analysis on the sequences shows that the selected primers are in the region of 2 non-synonymous mutations (S982A, A570D) and 2 synonymous mutations (C15279T, C16176T), with the best candidate being the one based on mutation A570D (Table 2).

We then further validate our results on 487 samples of other coronaviruses, obtained from the National Genomics Data Center (NGDC) [5]. All the candidate primers only target B.1.1.7 SARS-CoV-2, with no appearance in any other coronavirus sample. Further validation on 20,571 samples belonging to other non-corona viruses from the National Center for Biotechnology Information (NCBI) [START_REF] Stephen T Sherry | dbSNP: the NCBI database of genetic variation[END_REF], shows no appearance of the sequence in any other virus sample, providing further support for specificity of the obtained best candidate primers.

DISCUSSION

While the candidate primers identified by the proposed approach can be considered of high quality, it is worth it to compare them to other primer sets designed by conventional bioinformatics tools. A typical primer set would be represented by a sequence around a specific mutation, believed to be exclusive to the target virus variant.

In order to perform this comparison, we generate 21-bps sequences around mutations N501Y, A570D, D614G, P681H, T716I, S982A, D1118H, with 10 bps before and after the position of the mutation, using the canonical sequence of the B.1.1.7 variant: for example, mutation N501Y (A23063T) will correspond to sequence CCAACCCACT T ATGGTGTTGG. We then tested the frequency of appearance of these human-designed primers against the best primers found by the EA during the experimental evaluation. The results are reported in Table 3.

The 6 primer sets generated for the B.1.1.7 variant show almost the same frequency of appearance and similar specificity. Nevertheless, the sequence using mutation A570D is slightly better in frequency of appearance in the training set. The generated forward primers for the B.1.1.7 variant appear in 2,010 of 2,104 sequences, with an average frequency of 95.54%. A further analysis shows that only 2,014 of the sequences labeled as B.1.1.7 present 5 or more of the 7 studied mutations, which could point to an error in the annotation of the variant inside the GISAID dataset, or several extra mutations in the generated 21-bps sequences. If we consider as proper B.1.1.7 variants only sequences that show 5 or more of While the primers generated by our approach have a performance similar to primers obtained by conventional PCR design, it is extremely important to remark that this primers can only exist after the canonical sequence of a virus variant is defined. The canonical sequence represents the reference genome of a virus, and a considerable amount of effort from experts is necessary to reach a consensus on this sequence. Our approach, like other ML techniques for primer design [START_REF] Lopez-Rincon | Design of Specific Primer Set for Detection of B.1.1.7 SARS-CoV-2 Variant using Deep Learning[END_REF][START_REF] Lopez-Rincon | Classification and specific primer design for accurate detection of SARS-CoV-2 using deep learning[END_REF] has the advantage of automatically finding suitable primers in a matter of hours; but differently from similar ML approaches, it is much faster (about an order of magnitude faster on comparable hardware, from the experiments reported in the previously cited publications) and can handle a wider variety of constraints describing the desired features of a candidate primer.

As of February 5th, other SARS-CoV-2 variants of concern have been identified and are on the rise across the globe, such as the one originating from Brazil (P.1) [START_REF] Faria | Genomic characterisation of an emergent SARS-CoV-2 lineage in Manaus: preliminary findings[END_REF][START_REF]Emergence of SARS-CoV-2 B.1.1.7 Lineage -United States[END_REF] and the one generated in South Africa (B.1.351) [START_REF]SARS-CoV-2 Variants[END_REF][START_REF] Hadfield | Nextstrain: real-time tracking of pathogen evolution[END_REF][START_REF] Tegally | Emergence and rapid spread of a new severe acute respiratory syndrome-related coronavirus 2 (SARS-CoV-2) lineage with multiple spike mutations in South Africa[END_REF]. These new SARS-CoV-2 variants also carry the N501Y and D614G mutations, similarly to the B.1.1.7 variant. Thus, it is important to verify that the primers generated in this work are able to differentiate between the variants.

From the GISAID repository, we downloaded all available 326 sequences of the B.1.351 variant on January 7, 2021 and we downloaded all 28 non-repeated sequences of the P.1 variant on January 19, 2021. Next, we verified the frequency of appearance of the primers in samples from other variants. From Table 4, it is evident that while the primers found with our approach are exclusive to the B.1.1.7 variant, other sequences built around mutations (such as the one built around N501Y) are often also found in different variants. thus negatively impacting their specificity.

Finally, from a comparison between the frequency of appearance of reverse primers in UK samples (Table 5) and the frequency of appearance in other variants (Table 4), we are able to conclude that the best candidates for a primer set are sequences GCA CGT CTT GAC AAA GTT GAG as forward primer, based on mutation S982A, and GAG GTG CTG ACT GAG GGA AG as reverse primer. While the results are encouraging, it is important to remember that an in-silico validation is not enough to provide a final answer. The next step in the process will be to test the primers in the lab.

CONCLUSIONS

A wide variety of diagnostic tests have been used by high-throughput national testing systems around the world to monitor the SARS-CoV-2 infection [START_REF] Afzal | Molecular diagnostic technologies for COVID-19: Limitations and challenges[END_REF]. The arising prevalence of new SARS-CoV-2 variants such as B.1.1.7 has become of great concern, as most RT-PCR test to date will not be able to distinguish these new variants because they where not designed for such a purpose. Therefore, public health officials most rely on their current testing systems and their sequencing results to draw conclusions on the prevalence of new variants in their territories. An example of such cases has been seen in the UK, where the increase of the B.1.1.7 SARS-CoV-2 variant infection in their population was identified only through an increase in the S-gene target failure in their three target gene assay (N+, ORF1ab+, S-), coupled with sequencing of the virus and RT-PCR amplicons products [START_REF] Chand | Investigation of novel SARS-COV-2 variant: Variant of Concern 202012/01[END_REF]. Researchers believe that the Sgene target failure occurs due to the failure of one of the RT-PCR probes to bind as a result of the △69-70 deletion in the SARS-CoV-2 spike protein, present on B.1.1.7 [START_REF] Chand | Investigation of novel SARS-COV-2 variant: Variant of Concern 202012/01[END_REF]. This △69-70 deletion, which affects its N-terminal domain, has been recurrently occurring in different SARS-CoV-2 variants around the world [START_REF]Rapid increase of a SARS-CoV-2 variant with multiple spike protein mutations observed in the United Kingdom[END_REF][START_REF] Hadfield | Nextstrain: real-time tracking of pathogen evolution[END_REF][START_REF] Shu | GISAID: Global initiative on sharing all influenza data-from vision to reality[END_REF] and has been associated with other spike protein receptor binding domain changes [START_REF] Arambaut | Preliminary genomic characterisation of an emergent SARS-CoV-2 lineage in the UK defined by a novel set of spike mutations[END_REF]. Due to the likeliness of mutation in the S-gene, assays relying solely on its detection are not recommended, and a multiplex approach is required [START_REF] Afzal | Molecular diagnostic technologies for COVID-19: Limitations and challenges[END_REF][START_REF]Molecular assays to diagnose COVID-19: summary table of available protocols[END_REF][START_REF] Wang | Mutations on COVID-19 diagnostic targets[END_REF]. This is consistent with other existing designs, like CoV2R-3 in the S-gene [START_REF] Kim | The Progression of SARS Coronavirus 2 (SARS-CoV2): Mutation in the Receptor Binding Domain of Spike Gene[END_REF], that will also yield negative results for the B.1.1.7 variant, as the reverse primer sequence is in the region of mutation P681H. A more indepth analysis of S-dropout positive results can be found in Kidd et al. [START_REF] Kidd | S-variant SARS-CoV-2 is associated with significantly higher viral loads in samples tested by ThermoFisher TaqPath RT-QPCR[END_REF]. Given the concern for the increase in prevalence of the new SARS-CoV-2 B.1.1.7 variant and its possible clinical implication in the ongoing pandemic, diagnosing and monitoring the prevalence of such variant in the general population will be of critical importance to help fight its spread and develop new policies.

In this work, we propose an evolutionary approach to the design of primer sets. The approach is tested on the identification of the B.1.1.7 SARS-CoV-2 variant. The primer set developed using our approach proves to perform on par with primers found through other methods, and can be obtained in a reduced amount of time. We believe that the primer set can be used in a multiplexed approach in the initial diagnosis of COVID-19 patients, or used as a second step of diagnosis in cases already verified positive to SARS-CoV-2, to identify individuals carrying the B.1.1.7 variant. In this way, health authorities can better evaluate the medical outcome of this patients, and adapt or inform new policies that can help curve the rise of infections by variants of interest. Although the primer sets delivered by our automated methodology will still require laboratory testing to be validated, our approach can enable the timely, rapid, and low-cost operations needed for the design of new primer sets to accurately diagnose new emerging SARS-CoV-2 variants and other infectious diseases.

Figure 1 :

 1 Figure 1: PCR amplification of a section of interest in viral RNA. A primer uniquely identifies a section.

Figure 2 :

 2 Figure 2: Subsequence 𝐼 , representing a primer candidate, is uniquely identified by two integer values: index of the sample 𝑘 in the training set, and position inside its genome 𝑝.

21 𝑖=0𝑁

 21 (𝐼 (𝑖)) where 𝑁 (𝑏) = 1, if base 𝑏 is N 0, otherwise. (4) that takes into account the presence of N symbols in the sequence, indicating an error in the read. The ideal primer candidate should only contain A, C, G, or T values. The final term tackles the requirement of having a melting temperature 𝑇 𝑚 centered around 60 • . Specialized literature [34] provides an equation to compute 𝑇 𝑚 for a sequence 𝐼 : 𝑇 𝑚 (𝐼 ) = 81.5 + 16.6 * 𝑙𝑜𝑔 10 ([𝑁 𝑎+]) + 41 * C(𝐼 ) -600/𝑙 (𝐼 ) (5)

Table 1 :

 1 Candidate primers found during 20 runs of the proposed approach, and Primer3Plus output simulated on the canonical reference sequence 𝐸𝑃𝐼 _𝐼𝑆𝐿_601443.

	Sequence	Result
	GCACGTCTTGACAAAGTTGAG GAGGTGCTGACTGAGGGAAG
	TGGCAGAGACATTGATGACAC Left primer is unacceptable: High end self complementarity
	CAGAGACATTGATGACACTAC Left primer is unacceptable: Tm too low
	GGCAGAGACATTGATGACACT Left primer is unacceptable: High end self complementarity
	CCTCAAGGTATTGGGAACCTG CATCACAACCTGGAGCATTG
	TTGGCAGAGACATTGATGACA AGCAACAGGGACTTCTGTGC
	ACCTCAAGGTATTGGGAACCT CATCACAACCTGGAGCATTG
	CACACAACACATTTGTGTCTG Left primer is unacceptable: Tm too low/High end self complementarity
	TTCAGTGCATCGATATCGGTA	Left primer is unacceptable: High end self complementarity
	CTCAGACTAATTCTCATCGGC Left primer is unacceptable: Tm too low/High 3' stability
	TCAGACTAATTCTCATCGGCG Left primer is unacceptable: High 3' stability
	TCAGACTAATTCTCATCGGCG Left primer is unacceptable: High 3' stability
	GGCAGAGACATTGATGACACT Left primer is unacceptable: High end self complementarity
	CAGACTAATTCTCATCGGCGG Left primer is unacceptable: High 3' stability
	GTGATGTAGAAAACCCTCATC Left primer is unacceptable: Tm too low
	GTGATGTAGAAAACCCTCATC Left primer is unacceptable: Tm too low
	CTCAGACTAATTCTCATCGGC Left primer is unacceptable: Tm too low/High 3' stability
	CTCAGACTAATTCTCATCGGC Left primer is unacceptable: Tm too low/High 3' stability
	CTCATCTTATGGGTTGGGATT	GCCACACATGACCATTTCAC
	CCTTGCACGTCTTGACAAAGT GAGGTGCTGACTGAGGGAAG

Table 2 :

 2 Frequency of appearance in the training and test set, for each of the six candidate primers validated by Primer3Plus.

	Candidate Primer	Frequency in test set Frequency in training set Mutation
	GCACGTCTTGACAAAGTTGAG	0.9922	0.9893 T24506G (S982A)
	CCTCAAGGTATTGGGAACCTG	0.9922	0.9889 C16176T
	TTGGCAGAGACATTGATGACA	0.9922	0.9895 C23271A (A570D)
	ACCTCAAGGTATTGGGAACCT	0.9922	0.9888 C16176T
	CTCATCTTATGGGTTGGGATT	0.9910	0.9893 C15279T
	CCTTGCACGTCTTGACAAAGT	0.9922	0.9890 T24506G (S982A)

Table 3 :

 3 Frequency of appearance of the most significant mutations and the forward primer in the 8,293 sequences from the GISAID dataset.

	Sequence	B.1.1.7	Other Variants
		# samples (%) # samples (%)
	N501Y	1,985 (94.34%) 14 (0.02%)
	A570D	2,013 (95.67%) 1 (< 0.01%)
	D614G	2,096 (99.62%) 5,384 (78.96%)
	P681H	2,014 (95.72%) 1 (< 0.01%)
	T716I	2,005 (95.29%) 1 (< 0.01%)
	S982A	2,008 (95.43%) 0 (0%)
	D1118H	2,011 (95.57%) 0 (0%)
	GCACGTCTTGACAAAGTTGAG 2,011 (95.57%) 0 (0%)
	CCTCAAGGTATTGGGAACCTG 2,008 (95.43%) 0 (0%)
	TTGGCAGAGACATTGATGACA 2,014 (95.72%) 1 (< 0.01%)
	ACCTCAAGGTATTGGGAACCT 2,007 (95.38%) 0 (0%)
	CTCATCTTATGGGTTGGGATT	2,012 (95.62%) 1 (<0.01%)
	CCTTGCACGTCTTGACAAAGT 2,009 (95.48%) 0 (0%)
	Total Samples	2104	6137

Table 4 :

 4 Frequency of appearance of the characteristic mutations for the UK (B.1.1.7), South African (B.1.351), and Brazilian (P.1) variants in B.1.1.7 sequences (2,104), B.1.351 sequences (337), P.1 sequences (28) and sequences of other variants (6,808).

	Sequence	Other B.1.1.7 B.1351 P.1
	T1001I	0.01% 95.53% 0.00%	0.00%
	A1708D	0.00% 91.83% 0.00%	0.00%
	I2230T	0.01% 94.39% 0.00%	0.00%
	N501Y	0.04% 94.34% 99.70% 82.14%
	A570D	0.01% 95.67% 0.00%	0.00%
	P681H	0.01% 95.72% 0.30%	0.00%
	T716I	0.01% 95.29% 0.00%	0.00%
	S982A	0.00% 95.44% 0.00%	0.00%
	D1118H	0.00% 95.58% 0.00%	0.00%
	Q27stop	0.04% 90.64% 0.30%	0.00%
	R52I	0.00% 90.64% 0.00%	0.00%
	Y73C	0.00% 90.78% 0.00%	0.00%
	S235F	0.03% 95.58% 0.00%	0.00%
	GCACGTCTTGACAAAGTTGAG 0.00% 95.58% 0.00%	0.00%
	CCTCAAGGTATTGGGAACCTG 0.00% 95.44% 0.00%	0.00%
	TTGGCAGAGACATTGATGACA 0.01% 95.72% 0.00%	0.00%
	ACCTCAAGGTATTGGGAACCT 0.00% 95.39% 0.00%	0.00%
	CTCATCTTATGGGTTGGGATT	0.03% 95.63% 0.00%	0.00%
	CCTTGCACGTCTTGACAAAGT 0.00% 95.48% 0.00%	0.00%

Table 5 :

 5 Frequency of appearance of the forward and reverse primers found by the EA algorithm.

	Forward Primer	Frequency Reverse Primer	Frequency Average
	GCACGTCTTGACAAAGTTGAG 95.58%	GAGGTGCTGACTGAGGGAAG 99.71%	97.65%
	CCTCAAGGTATTGGGAACCTG 95.44%	CATCACAACCTGGAGCATTG 98.62%	97.03%
	TTGGCAGAGACATTGATGACA 95.72%	AGCAACAGGGACTTCTGTGC 99.62%	97.67%
	ACCTCAAGGTATTGGGAACCT 95.39%	CATCACAACCTGGAGCATTG 98.62%	97.01%
	CTCATCTTATGGGTTGGGATT	95.63%	GCCACACATGACCATTTCAC 99.81%	97.72%
	CCTTGCACGTCTTGACAAAGT 95.48%	GAGGTGCTGACTGAGGGAAG 99.71%	97.60%

https://www.nature.com/scitable/definition/primer-305/

The repository address will be disclosed after review, to prevent identification of the authors, as per double-blind peer-review rules.


Alejandro Lopez-Rincon
email: a.lopezrincon@uu.nl

Daphne S Roozendaal

Hilde M Spierenburg

Asta L Holm

Renee Metcalf

Paula Perez-Pardo

Aletta D Kraneveld

Alberto Tonda
email: alberto.tonda@inrae.fr

Modelling Asthma Patients' Responsiveness to Treatment Using Feature Selection and Evolutionary Computation

Keywords: Omalizumab, Machine learning, Evolutionary computation, Asthma, Mathematical model

come

Introduction

When patients are treated with medical drugs, it is possible to observe a variation in their gene expression. Such transcriptional variation can potentially be correlated with the responsiveness to treatment, and this relationship can be described through a system of Ordinary Differential Equations (ODE).

Such a model would be important not only to explain the differences in treatment outcome, but also to provide indications to medical personnel on how to improve the therapy. For example, if a gene is shown to be overexpressed in non-responsive patients, with respect to responsive patients, doctors might devise a new therapy, combining the current treatment with substances that lower the expression of that particular gene.

Nevertheless, obtaining such an ODE system is not a straightforward process. The genes potentially correlated to drug response number in the tens of thousands, with possible complex interconnections in expression levels. Not only that, but even when a limited number of genes is identified, finding satisfying values of the parameters for the ODE system describing their interactions is a complex optimization problem that cannot be tackled through gradient-based techniques.

In this paper, we propose a novel methodology to obtain compact mathematical models describing the correlation between gene expression levels and responsiveness to treatment. The methodology combines a technique for feature selection [START_REF] Lopez-Rincon | Automatic discovery of 100-mirna signature for cancer classification using ensemble feature selection[END_REF][START_REF] Lopez-Rincon | Machine learning-based ensemble recursive feature selection of circulating mirnas for cancer tumor classification[END_REF], able to identify a small set of genes highly correlated with treatment outcome, and a state-of-the-art evolutionary optimizer [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] in order to find good values for the ODE system characterizing their interaction.

The presented approach is tested on real-world data from N=40 patients affected by moderate-to-severe asthma, treated with the recent anti-IgE drug Omalizumab (30 responsive and 10 non-responsive). The results show that the methodology is effective in identifying 5 genes that are highly correlated with responsiveness to treatment, and it is able to deliver an ODE system that can reliably describe their interaction, explaining the responsiveness of patients to the Omalizumab treatment.

Background

In this section, we introduce the minimal notions that are necessary to introduce the scope of our work.

Feature selection

In machine learning (ML), feature selection (FS) is defined as the process of identifying the features of a data set in order to obtain a minimal, informative subset. Features may not be part of this subset for two main reasons: they might be unrelated to the underlying nature of the problem, just adding noise; they might be heavily correlated with other features, adding no relevant information for the task. Applications range from face recognition [START_REF] Vignolo | Feature selection for face recognition based on multi-objective evolutionary wrappers[END_REF] to medicine [START_REF] Zhou | Multi-objective based radiomic feature selection for lesion malignancy classification[END_REF], and approaches can be divided into two categories [START_REF] Guyon | An introduction to variable and feature selection[END_REF]: filters that score features according to a criterion (often a statistical test); and recursive procedures (forward or backwards) that attempt to reduce the features to a small set of non-redundant ones [START_REF] Lewis | The characteristic selection problem in recognition systems[END_REF][START_REF] Chien | On the generalized karhunen-loève expansion (corresp.)[END_REF].

In the scope of this work, we focus on recursive FS algorithms, in particular Recursive Ensemble Feature Selection (REFS). The method is a variation of Recursive Feature Elimination (RFE) [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF] that scores the features in a 10-fold cross-validation scheme, using 8 different classifiers: gradient boosting, passive aggressive classifier, logistic regression, Support Vector Machine classifier (SVC), random forest, Stochastic Gradient Descent (SGD), ridge classifier and bagging. The lowest scoring features are removed from the analysis and the process is repeated until the overall classification accuracy drops below a given threshold. The use of an ensemble of classifiers reduces the effects of the inherent bias in each ML algorithm, thus delivering a more objective feature ranking. This technique has been applied successfully for problems involving both mRNA [START_REF] Lopez-Rincon | Machine learning-based ensemble recursive feature selection of circulating mirnas for cancer tumor classification[END_REF] and miRNA [START_REF] Lopez-Rincon | Automatic discovery of 100-mirna signature for cancer classification using ensemble feature selection[END_REF], featuring number of variables ranging from 1,046 to 54,675.

Omalizumab treatment for asthma patients

Omalizumab is the first humanized monoclonal antibody that binds to immunoglobuline E (anti-IgE) prescribed to patients with moderate-to-severe allergic asthma who do not respond to inhaled corticosteroids and long-acting β2agonist bronchodilators. Omalizumab works by specifically binding free serum IgE [START_REF] Thomson | Omalizumab: clinical use for the management of asthma[END_REF], which characterizes allergic asthma. IgE binds to high affinity receptors (Fc RI) expressed on effector cells such as basophils and mast cells, but also on other immune cells like eosinophils, thereby triggering an inflammatory cascade through the release of inflammatory mediators [START_REF] Hamelmann | The rationale for treating allergic asthma with anti-ige[END_REF]. By binding to the Fc region of IgE and forming IgE-antibody complexes, Omalizumab prevents the binding of allergen specific IgE to Fc RI and the subsequent inflammatory allergy reaction [START_REF] Thomson | Omalizumab: clinical use for the management of asthma[END_REF]. Consequently, Fc RI expression is reduced, leading to less immune activation. Furthermore, Omalizumab decreases eosinophil numbers found in the airway of asthmatic patients, although the mechanism through which this effect is achieved is not entirely understood [START_REF] Thomson | Omalizumab: clinical use for the management of asthma[END_REF]. Eosinophils and their derived proinflammatory mediators are major contributors to airway inflammation and damage [START_REF] Possa | Eosinophilic inflammation in allergic asthma[END_REF]. Omalizumab's ability to combat long-term airway remodeling is still under investigation [START_REF] Hamelmann | The rationale for treating allergic asthma with anti-ige[END_REF].

It is important to understand the mechanism of action of Omalizumab treatment and to research the differences in responsiveness. Potential transcriptional variations between responders (R) and non-responders (NR) to Omalizumab can function as predictive biomarkers in the future. A recent study by [START_REF] Upchurch | Whole blood transcriptional variations between responders and non-responders in asthma patients receiving omalizumab[END_REF], with accession number GSE134544 at gene expression omnibus (GEO) [START_REF] Upchurch | Whole blood transcriptional variations between responders and non-responders in asthma patients receiving omalizumab[END_REF] investigated whole blood transcriptomes of moderate-to-severe asthma patients (N = 40; 30 responders (R) and 10 non-responders (NR)), over the course of Omalizumab treatment. Blood was collected at day 0, 7, 42, 98, 182 where the treatment started at day 7, and day 0 marks one week before the treatment.

Total RNA was isolated from whole blood, and all samples passing quality control were then amplified, (biotin-)labelled and hybridized to Illumina HT-12 V4 BeadChips (Illumina). Subsequent differential gene expression analysis was performed using Welch's T-test for comparisons between R and NR, and transcriptional changes within each group were assessed using a paired T-test. However, direct comparisons between R and NR did not provide sufficient gene lists after multiple testing corrections. Therefore, whole blood mRNA signature differences between groups were characterized using a gene cluster strategy, or modular-level type analysis adopted from studies by Chaussabel et al. [START_REF] Chaussabel | A modular analysis framework for blood genomics studies: application to systemic lupus erythematosus[END_REF] and Banchereau et al. [START_REF] Banchereau | Personalized immunomonitoring uncovers molecular networks that stratify lupus patients[END_REF]. Hierarchical clustering of genes revealed 8 similarly expressed transcript clusters in R and NR (i.e. protein synthesis (1); T cell/NK cell/ cytotoxicity (2); hematopoiesis (3); cell cycle control/proliferation (4); T cell regulation and activation (5); monocytes [START_REF] Chairakaki | Plasmacytoid dendritic cells drive acute asthma exacerbations[END_REF]; glucose metabolism [START_REF] Chaussabel | A modular analysis framework for blood genomics studies: application to systemic lupus erythematosus[END_REF] and inflammation ( 8)). Of these, cluster 2 and 7 in R were reported to be higher, while clusters 3 and 8 were higher in NR, suggesting that clusters 2, 3, 7 and 8 can be used as predictors of response to Omalizumab treatment. These clusters, combined, contain a total of 1,776 genes.

Due to the unpredictability of an asthma patient's responsiveness to Omalizumab, there have been multiple studies into finding a reliable biomarker that can act as a predictor. A recent study [START_REF] Huang | Endotypes of severe allergic asthma patients who clinically benefit from anti-ige therapy[END_REF] reported that interleukin (IL)-9, IL-13, IL-25, IL-33 and thymic stromal lymphopoietin (TSLP) levels were significantly higher in R compared to NR. Most R were of a high type-2 cytokine endotype whereas only one NR was [START_REF] Huang | Endotypes of severe allergic asthma patients who clinically benefit from anti-ige therapy[END_REF]. Data from the INNOVATE trial of Omalizumab for severe persistent asthma [START_REF] Humbert | Benefits of omalizumab as add-on therapy in patients with severe persistent asthma who are inadequately controlled despite best available therapy (gina 2002 step 4 treatment): Innovate[END_REF] was analyzed by [START_REF] Bousquet | Predicting and evaluating response to omalizumab in patients with severe allergic asthma[END_REF], and this study found that of the biomarkers recorded by INNOVATE, only baseline total IgE levels were a predictor of efficacy. However, pooled analysis showed that treatment was effective in some cases irrespective of IgE levels. It was concluded that the most meaningful measure of responsiveness was the physician's overall assessment [START_REF] Bousquet | Predicting and evaluating response to omalizumab in patients with severe allergic asthma[END_REF]. An extensive look into gene expression relating to asthma patients vs controls in different tissue types, disease severity and response to allergens and corticosteroid treatment in several datasets reported multiple gene signatures and pathways [START_REF] Alrashoudi | Gene expression data analysis identifies multiple deregulated pathways in patients with asthma[END_REF]. Although this may prove useful in explaining Omalizumab response, it revealed no significant gene overlap.

Proposed Approach

We present a new approach to obtain compact, human-readable mathematical models to explain responsiveness to treatment in patients. The methodology first applies feature selection to identify a small set of relevant genes, and then uses state-of-the-art evolutionary optimization to find the parameters of an ODE system that describes the relationship between gene expression levels and patient's responsiveness to treatment.

Feature Selection

In a first step, our objective is to select the most meaningful genes to correctly predict and model patients' responsiveness to treatment. We apply the REFS algorithm, which uses the feedback of an ensemble of classifiers to rank each feature depending on its usefulness for the process of classification. Then, the lowest-scoring features are removed, and the classification/ranking is repeated, until the average classification accuracy falls below a user-defined threshold.

Mathematical Modeling

Once a small subset of the genes is identified, we create a mathematical model that interconnects the gene expression given the values at different time points. While other solutions to model the correlation between the gene expression values and responsiveness are possible, like black-box machine learning, white-box models are preferred by practitioners, as they are commonly considered more interpretable. We assume to have whole blood mRNA expression at different points in time available, and we will consider the average value for R and NR over all samples at each point in time, for each category of patients. Thus, the mathematical model we propose is a system of ordinary differential equations (ODEs) in the form of Eq. 1:

dg 0 dt = -k 0 u(t) + α 0 e -β0t , ..., dg n dt = -k n u(t) + α n e -βnt , u(t) = K g0 g 0 + ... + K gn g n , (1)
where g 0 to g n will be the most important genes, k i , α i , β i and K gi are coefficients calculated by each gene and u(t) is an unknown function that interconnects the gene expression of all the genes, as to consider a relationship between all the variables.

Then, to solve the model, we use Euler's numerical method which transforms system 1 into:

g t 0 = g t-1 0 + t dg 0 dt , ..., g t n = g t-1 n + t dg n dt , u(t) = K g0 g t 0 + ... + K gn g t n . (2)
u(t) models the interconnection between different gene expressions, that we hypothesize exists to avoid trivial assumptions of independence. As the problem is not treatable resorting to classical gradient-based techniques, it is necessary to use state-of-the-art stochastic optimization, such as CMA-ES [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF], to find satisfying values for k i , α i , β i and K gi .

Given the measurement of gene expressions at different instants t = {t 0 , ..., t N }, from Eq 2 we can define the cost function to be minimized by CMA-ES as:

error = n i=0 |g t CRi -g t Ri | + |g t CN Ri -g t N Ri | (3)
where g t CRi is the average gene i expression calculated at time t for R, g t Ri is the average gene i expression measured t time t for responders, g t CN Ri is the average gene i expression calculated at time t for non-responders and g t N Ri is the average gene i expression measured t time t for non-responders.

Experimental Evaluation

All the necessary code for the experiments has been developed in Python, using the scikit-learn package [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF] for machine learning, and the cma package for CMA-ES. The code is hosted on the open GitHub repository: 3

Data

Moderate-to-severe asthma patients were prescribed Omalizumab, based on the manufacturer's dosing table. Patient blood was collected at day 0, 7, 42, 98, 182 where the treatment started at day 7 with Omalizumab in 40 patients: 30 R and 10 NR. For each patient, for each sample, for each instant of time, the dataset contains information about 28,402 gene expression levels. All data was used as provided, at GEO accession code GSE134544 [START_REF] Upchurch | Whole blood transcriptional variations between responders and non-responders in asthma patients receiving omalizumab[END_REF].

Feature Selection

Running the REFS algorithm previously described 10 times, we identified a set of 5 features (out of [START_REF] Thomas | Ciliary dysfunction and ultrastructural abnormalities are features of severe asthma[END_REF]402). This compact set can predict the Omalizumab responsiveness in patients with a mean accuracy of 0.975 in a 10-fold cross-validation, considering the binary classification problem (R/NR) with all classifiers in the REFS ensemble. As the REFS process is stochastic, it was iterated 10 times and the feature set corresponding to the highest peak in accuracy was selected (see Fig. 1). From the figure, it is interesting to notice how using all 28,402 features actually provides a lower mean classification accuracy (0.703). Classification algorithms, usually exploiting optimization heuristics, often show a lower performance when asked to explore a larger feature search space.

The resulting most significant features uncovered by the presented algorithms are ILM N 3286286; ILM N 1775520; ILM N 1656849; ILM N 1781198 and ILM N 1665457 (Fig. 2). Details of the Illumina probes are further specified in Table 1. To further validate the selected features, we computed the area under the curve (AUC) and receiver operating characteristic (ROC) curve in a 10-fold cross-validation, using just the selected features, testing all classifiers in the REFS ensemble. The best AUC, 0.99, was obtained using Passive Aggressive classifier. This result is considered as an excellent diagnostic accuracy (AUC 0.9-1.0) by specialists of the field [START_REF] Šimundić | Measures of diagnostic accuracy: basic definitions[END_REF][START_REF] Mandrekar | Receiver operating characteristic curve in diagnostic test assessment[END_REF].

Mathematical Modeling

Considering we find 5 gene expressions to be the most meaningful, and using the template defined in Eq. 1, we can write the specific model linking gene expression levels to responsiveness to Omalizumab as:

dg 0 dt = -k 0 u(t) + α 0 e -β0t , dg 1 dt = -k 1 u(t) + α 1 e -β1t , dg 2 dt = -k 2 u(t) + α 2 e -β2t , dg 3 dt = -k 3 u(t) + α 3 e -β3t , dg 4 dt = -k 4 u(t) + α 4 e -β4t , u(t) = K g0 g 0 + K g1 g 1 + K g2 g 2 + K g3 g 3 + K g4 g 4 , (4)
where we will need to find 20 parameters. We run CMA-ES with a σ = 0.01 and λ = 1000 with default stop conditions. Using CMA-ES with the cost function in Eq. 3, with t = 0.25, we find the following values; k 0 = 0.5313, k 1 = 0.6440, k 2 = 0.4899, k 3 = 0.5504, k 4 = -0.4702, α 0 = -0.0001, β 0 = -0.0168, α 1 = -0.0386, β 1 = 0.9277, α 2 = -0.0071, β 2 = 0.0733, α 3 = -0.0176, β 3 = 0.1422, α 4 = 0.0029, β 4 = 0.0079, K g0 = 0.2598, K g1 = 0.1078, K g2 = 0.8069, K g3 = -0.5708, K g4 = 0.0565 with an error = 1.6036 as the best of 20 runs.

As a baseline comparison, we also tested the scipy [START_REF] Virtanen | SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF] implementation of the Nelder-Mead optimization algorithm [START_REF] Gao | Implementing the nelder-mead simplex algorithm with adaptive parameters[END_REF], currently considered among the state-of-the-art for gradient-free optimization, for 20 runs. The best run of the Nelder-Mead algorithm yields a solution with error = 4.3143, of lower quality than that of CMA-ES.

From the results in Fig. 4, we can see that the model approximates considerably the behaviour of the genes, given the parameters and the initial values only; with a clear exception of LOC100131780. Although the increase of responders and decrease of non-responders functions are reflected, the amplitude does not match, therefore suggesting the necessity of increasing the degree of the answer or a more precise ODE solver.

Finally, from function u(t) reported in Fig. 5, it is possible to notice that, just from its initial value u(0), it is already possible to differentiate responders from non-responders, predicting the outcome of the treatment before its beginning.

Discussion and Conclusions

Severe-to-moderate asthma patients appear to respond differently to the biological anti-IgE treatment Omalizumab. Total mRNA sequencing of whole blood from R and NR of Omalizumab treatment can help appreciate the differences at a transcriptional level. In this study, the GSE134544 dataset (that was recently published by Upchurch et al. [START_REF] Upchurch | Whole blood transcriptional variations between responders and non-responders in asthma patients receiving omalizumab[END_REF]), was analyzed with the use of machine learning-based REFS. This novel study reveals 5 genes that are highly relevant in predicting Omalizumab responsiveness in asthma patients. In addition, we created a mathematical model to approximate the interdependence of the most significant genes to explain the effect of Omalizumab treatment using evolutionary computation.

Investigating gene function corresponding to the 5 mRNA (Table 1) predictors found in this study may illuminate new pathways involved in allergic asthma and the mechanism of Omalizumab resistance. For example, CCDC113 is vital for ciliogenesis and when knocked down, causes a reduction in cilium formation [START_REF] Firat-Karalar | Proteomic analysis of mammalian sperm cells identifies new components of the centrosome[END_REF] which is previously related to severe asthma [START_REF] Thomas | Ciliary dysfunction and ultrastructural abnormalities are features of severe asthma[END_REF]. Also, this gene was found to be overexpressed in asthma patients as compared to controls [START_REF] Pandey | A nasal brush-based classifier of asthma identified by machine learning analysis of nasal rna sequence data[END_REF], and in NR as compared to R.

CLEC4C is a marker for plasmacytoid dendritic cells (pDCs) subtypes [START_REF] Murray | Clec4c gene expression can be used to quantify circulating plasmacytoid dendritic cells[END_REF][START_REF] Murray | Risks for cold frequency vary by sex: role of asthma, age, tlr7 and leukocyte subsets[END_REF], which both have been implicated in driving acute asthma exacerbations [START_REF] Chairakaki | Plasmacytoid dendritic cells drive acute asthma exacerbations[END_REF][START_REF] Vroman | Dendritic cell subsets in asthma: impaired tolerance or exaggerated inflammation?[END_REF] and have shown to have a tolerogenic effect in asthma by inducing Treg cell differentiation [START_REF] Vroman | Dendritic cell subsets in asthma: impaired tolerance or exaggerated inflammation?[END_REF]. CLEC4C is overexpressed in R when comparing against NR. PPP1R3D, a gene that codes for a subunit of PP1 (protein phosphatase 1) [START_REF] Armstrong | Ppp1r6, a novel member of the family of glycogen-targetting subunits of protein phosphatase 1[END_REF] regulates protein serine/threonine phosphatase activity and was found to be a causal key driver for acute peanut allergic response [START_REF] Watson | Integrative transcriptomic analysis reveals key drivers of acute peanut allergic reactions[END_REF]. Genes coding for other subunits of PP1 have been shown to be upregulated in asthmatic patients (PPP1R16A) or, more specifically, in corticosteroid resistant patients (PPP1R15A) [START_REF] Alrashoudi | Gene expression data analysis identifies multiple deregulated pathways in patients with asthma[END_REF]. This gene is underexpressed in R, when compared against NR.

More unexpected in its relation to Omalizumab responsiveness in severe-tomoderate asthma is SLC26A8, which is most commonly associated with sperm motility and mutations that can cause male infertility [START_REF] Dirami | Missense mutations in slc26a8, encoding a sperm-specific activator of cftr, are associated with human asthenozoospermia[END_REF]. It was also found to be upregulated in patients with severe asthma as compared to healthy controls [START_REF] Bigler | A severe asthma disease signature from gene expression profiling of peripheral blood from u-biopred cohorts[END_REF]. Lastly, LOC100131780 is proprietary of Illumina and we, unfortunately, could not find a perfect gene match for this probe. The sequence did, however, overlap partially with DNAI1, which codes for dynein axonemal intermediate chain 1. The gene is strongly linked to primary ciliary dyskinesia (DNAI1 -Dynein Intermediate Chain 1, Axonemal -Homo Sapiens (Human) -DNAI1 Gene & Protein), which can cause respiratory infections and breathing problems (NHLBI).

An attempt to understand transcriptional variations between responders and non-responders on Omalizumab treatment was made by Upchurch et al. 2020, who recently analyzed changes in 8 gene clusters. In Omalizumab responders, the T cell/natural killer (NK) cell/cytotoxicity gene cluster and the glucose metabolism gene cluster were higher, whereas gene clusters involved in hematopoiesis and inflammation appeared to be higher in non-responders. 19 genes of the 8th (inflammatory) gene cluster were specifically annotated in their paper [START_REF] Upchurch | Whole blood transcriptional variations between responders and non-responders in asthma patients receiving omalizumab[END_REF]; however, none of these genes corresponded with the 5 most relevant genes we found with our methodology.

By setting a goal of predicting responsiveness instead of investigating transcriptional changes, the amount of data and complexity are reduced, since only the pre-testing time point needs to be considered. Using REFS provides the benefit of illuminating possible gene interactions without bias, as opposed to using a single algorithm or performing significance tests and clustering up-or downregulated genes on the basis of known pathway function. The overall benefit of machine learning over basic statistics is that it is able to find predictive patterns which do not rely on assumptions about the data-generating system. Furthermore, it avoids the problem of eliminating significant mRNAs by not needing to take into account false discovery rate (FDR correction). Therefore, the technique used in this paper greatly reduced the number of meaningful mRNAs (5) compared to Upchurch et al. 2020 [START_REF] Alrashoudi | Gene expression data analysis identifies multiple deregulated pathways in patients with asthma[END_REF]776).

Further clinical testing and a replication dataset could reveal whether the 5 mRNAs proposed in this paper can reliably predict moderate-to-severe asthma patients' response to Omalizumab treatment on a larger scale. Implementing this sort of pre-treatment testing can both reduce the cost of asthma treatment, as Omalizumab is a relatively expensive drug, and prevent unnecessary and unproductive treatment time. Although currently, we limit our study to Omalizumab, it could be applied to similar studies.

This study not only proposes 5 specific genes that are transcriptionally relevant to predict Omalizumab responsiveness in moderate-to-severe asthma patients, but also puts forward a novel technique that aims to reduce the necessary information to the smallest set of whole blood mRNAs. The presented methodology elucidates the power of machine learning versus more general univariate/multivariate statistical analysis strategies.

Fig. 1 .

 1 Fig. 1. The results of 10 runs of the REFS algorithm for the classification of Omalizumab responsiveness in allergic asthma patients. The x axis cuts at 5 variables, in correspondence with the highest peak.

Fig. 2 .

 2 Fig.2. Heatmap for the normalized gene expression for the 5 selected genes in all of the samples. Interestingly, samples of the two classes can be visually separated by just looking at the different normalized gene expression of the samples in the two groups.

Fig. 3 .

 3 Fig. 3. ROC Curve in a 10-fold cross validation using Passive Aggressive classifier (the most effective in the REFS ensemble for this particular problem) for the 5 selected genes.

Fig. 4 .

 4 Fig. 4. Genes expression modeling and measured values for the 5 selected features at t = 0, 7, 42, 98, 182.

Fig. 5 .

 5 Fig. 5. Complementary signal u(t) = Kg 0 g0 + Kg 1 g1 + Kg 2 g2 + Kg 3 g3 + Kg 4 g4 that interconnects the expression values of all of the genes.

Table 1 .

 1 Information on the selected features to predict the responsiveness of the Omalizumab treatment.

	Illumina probe Corresponding gene Sequence
	ILMN 3286286 LOC100131780	GAGATTGCGAAACTGGACAAACTGC TGAACCTGGACAGGGGCCAGGGCTG
	ILMN 1775520 CCDC113	GGACATGAGAACATATTTCCAAGACA GAGGATTCTATGGGGACGGGTCAC
	ILMN 1656849 SLC26A8	TGGGCGTATTGGGTTTGGGCTTCATTG CCACTTACCTTCCGGAGTCTGCA
	ILMN 1781198 PPP1R3D	GGCCTTCACTGCTACGCCCCTGGCCCC AAAACAGAGAGCAAGACAGTTGT
	ILMN 1665457 CLEC4C	GTGGTTCCAGTTGAAGGTCTGGTCCAT GGCAGTCGTATCCATCTTGCTCC

https://github.com/steppenwolf0/modelingEvolutionaryComputation