Abstract
We investigate the optimal vaccination and screening strategies to minimize human papillomavirus (HPV) associated morbidity and the interventions cost. We propose a two-sex compartmental model of HPV-infection with time-dependent controls (vaccination of adolescents, adults, and screening) which can act simultaneously. We formulate optimal control problems complementing our model with two different objective functionals. The first functional corresponds to the protection of the vulnerable group and the control problem consists of minimizing the cumulative level of infected females over a fixed time interval. The second functional aims to eliminate the infection, and, thus, the control problem consists of minimizing the total prevalence at the end of the time interval. We prove the existence of solutions for the control problems, characterize the optimal controls, and carry out numerical simulations using various initial conditions. The results and properties and drawbacks of the model are discussed.
1. Introduction
Human papillomavirus (HPV) is the leading etiological factor for the development of cervical cancer. According to the World Health Organization (WHO) report, cervical cancer is the fourth most frequent cancer in women with approximately 570,000 new cases reported in 2018 [1]. This fact makes the design, implementation, and maintenance of effective preventing policies against cervical cancer a highly relevant public health problem. This consideration motivates research on the design and development of cost-effective vaccination and screening programs against HPV infection.
On October 5th of 2018, the US Food and Drug Administration (FDA) approved the nonavalent vaccine Gardasil-9 against HPV infection for use in women and men aged 27 to 45 [2]. Earlier, the FDA approved the vaccine for persons between the ages of 9 to 26. Gardasil-9 protects against HPV types 6, 11, 16, and 18 (that are already targeted by the quadrivalent HPV vaccine Gardasil) and against the next five most common oncogenic viral types, namely, HPV 31, 33, 45, 52, and 58. These nine HPV types are responsible for the majority of HPV-associated diseases [3]. In addition to the nonavalent and quadrivalent vaccines, there is also an FDA-approved bivalent vaccine (Cervarix) that targets HPV types 16 and 18.
The introduction of HPV vaccines is changing the epidemiology of cervical cancer and other HPV-associated diseases each year. Since HPV-associated morbidity and mortality are usually considerably higher in females than in males, in the majority of countries the primary target group of HPV vaccination programs is adolescent school girls aged 9–14 [4]. However, considering that Gardasil-9 has been licensed for a broader age range in both males and females, it is essential to analyze the potential health benefits of extending vaccination to males, as well as older individuals, versus a higher economic cost of such programs.
To examine the dynamics of HPV spread and control, numerous mathematical models have been proposed [5–15]. While these and other recent studies analyzed different strategies for HPV control and brought important insight into the problem, the majority of them did not take into account the optimality of these interventions. The study of optimality, however, is of particular relevance in the allocation of limited public health resources. In this context, the application of optimal control theory has proven to be an important tool for evaluating and optimizing various detection, prevention, therapy, vaccination, and other intervention programs [16].
So far, to the best of authors’ knowledge, for the case of HPV infection only a few studies have been done using the optimal control framework to assess the impact of vaccination programs [17, 18]. Brown and White [17] explored how to target vaccination in the UK. This study highlights the importance of including a catch-up vaccination policy in order to control the spread of the infection. In a recently published paper, Malik et al. [18] analyzed the optimal control strategies for a vaccination program administering the bivalent, quadrivalent, and nonavalent HPV vaccines to the female population. In particular, they explored scenarios where the three vaccines are used simultaneously compared with the case where the bivalent and quadrivalent vaccines were initially used and then, during the program, one or both of them are replaced by the nonavalent vaccine.
Since there is evidence that switching from the quadrivalent to the nonavalent is likely to be cost-saving [3], in this work we focus on the application of the nonavalent vaccine. The aim is to use optimal control theory to gain insight into the best combination of vaccination and screening to reduce the spread of HPV infection, as well as the cost of the intervention strategy. To address this issue, we construct a mathematical model that allows vaccination prior to and after sexual initiation in both males and females. Moreover, we also include into the model a screening program that allows detection of HPV-associated diseases.
In general, public health interventions have two basic objectives: protection of population groups that are at considerably higher risk of contracting or developing severe disease and total elimination of the infection when it is potentially possible. Accordingly, in this paper, we consider both of these objectives formulating two optimal control problems. Our ultimate goal is to explore where control policies aimed at the protection of the vulnerable group and at the eradication of the infection differ.
The organization of this paper is as follows: The next section is devoted to formulation and analysis of an epidemic model for the transmission dynamics of HPV infection in a heterosexual population. In Section 3, we propose optimal control problems incorporating two objective functionals and simultaneous application of several controls into our model. In Section 4, we prove the existence of solutions and characterization of the optimal controls via the Pontryagin maximum principle. Section 5 contains the results of numerical simulations of the model and the profiles of the optimal controls. The last section contains the conclusions and a discussion of the obtained results.
2. Model Formulation and Analysis
In this study, we propose and study a compartmental model for the transmission dynamics of the HPV types targeted by the nonavalent vaccine Gardasil-9. The model classifies the hosts’ population at time denoted by according to gender and infection status. We subdivide the total female population into four mutually exclusive compartments: unvaccinated susceptibles , vaccinated susceptibles , and infectious females unaware and aware of their infection; thusThe total male population consists of three compartments: unvaccinated susceptibles , vaccinated susceptibles , and infectious males , so that
We assume that individuals enter the sexually active population as singles at a constant rate with a sex ratio of 50:50. Females and males leave the population by death or ceasing sexual activity at per capita rates and , respectively. In this study, we distinguish between vaccinations prior and after sexual initiation. In particular, a fraction of females and a fraction of males are vaccinated before they enter the sexually active class and thus are recruited into their vaccinated compartment. Moreover, the susceptible sexually active females and males are vaccinated at rates and per unit of time, respectively. The vaccine reduces the force of infection by a factor ; thus, the vaccine is completely effective when and ineffective when . Therefore is the vaccine effectiveness. The vaccine-induced immunity wanes at a rate . Both males and females can clear the infection naturally with per capita rates and , respectively. No permanent immunity is assumed after that.
Susceptible females become infected at a rate following effective contact with an infectious male. After the acquisition of HPV infection, a fraction of the infected females may develop symptoms and become aware of their infection entering the class. The remaining infected females do not realize their infection and move to the unaware class . However, the screening program may allow the unaware infected females to detect their infection at a rate .
Susceptible males can be infected by unaware infected females at a rate and by aware infected females at a rate . We assume that since a female conscious of her infection can take precautions like condom use to reduce the possibility of transmission. Since HPV infection is transmitted (in most cases) sexually, we model the transmission via the standard incidence [22].
The model resulting from these assumptions is governed by the following system of differential equations:where all the parameters are assumed to be nonnegative.
We are interested in the dynamics of model (3) over a finite time interval . Therefore, we assume that the population for both sexes is constant with values and . This assumption is reasonable since the epidemic process is considerably faster than the demographic one over a short period of time.
Under the constant population assumption, we can express system (3) in terms of the proportionsOmitting the bars above the variables for convenience, our normalized system becomes
Please note that system (5) always has a unique disease-free equilibrium given by
In order to compute the basic reproduction number, we need the following components. During his infection period, , an infectious male produces on average infections in females. Therefore, is the infection transfer from infectious males to females. An infectious female aware of her infection produces on average infections during her infectious period . As a result, is the infection transfer from the aware infectious females to males. Analogously, is the infection transfer from the unaware infectious females to males. Moreover, a fraction of unaware infected females becomes aware of the infection via screening; thus is the infection transfer from females to males.
Using the next-generation matrixit can be shown that the basic reproduction number is the geometric mean of the terms and ; that is,As a consequence of heorem 2 in [23], the disease-free equilibrium is locally asymptotically stable if and unstable if .
3. The Optimal Control Problems
Five controls are possible in this modeling framework; hence, we formulate optimal control problems incorporating these five controls, namely:(i) and are the controls representing the fraction of females and males, respectively, that are vaccinated prior to sexual initiation(ii) and are the controls representing the vaccination rates of susceptible sexually active females and males, respectively(iii) is the control representing the screening rate for females
Because resources of a health system are limited, we assume that all the controls are subject to constraintsWith these controls and taking into consideration the constant population size assumption, we can write model (5) as the following system of five equations:complemented by initial conditionsThe set of admissible controls is the set of Lebesgue measurable functions which for almost all satisfy constraints (13).
We seek the optimal controls that reduce the incidence and prevalence of the infection and also the cost of using a public health intervention over years. In the case where the elimination of HPV infection is not possible due to social or economic constraints, the public health system can focus, instead, on the protection of the vulnerable group, that is, on minimizing the number of cervical cancer cases and related diseases in females over the time interval . For this case, the objective functional is Secondly, having the aim of the ultimate elimination of the infection in mind, we consider the problem of minimizing the total prevalence at the terminal time moment . In this context, the objective functional takes the following form:In these functionals, the parameter balances the probability of developing cervical cancer and its treatment cost, whereas reflects the cost of the long-term presence of HPV infection in the population. The meaning of the parameters , , and in both functionals is the same. Specifically, the weight parameter represents the cost associated with the vaccination of school girls and boys prior to their sexual initiation. The weight parameter represents the cost of administering the vaccine in sexually active individuals. Finally, represents the cost of screening and follow-up after an abnormal Pap test.
We are assuming that the cost of applying health policies increases as a nonlinear function of the controls. The quadratic terms penalize high levels of control administration in comparison with the cost of low levels. Due to its mathematical convenience, the sum of the weighted squares of the controls is probably the most common form of the objective functional in the literature: for functionals of this type, it is possible to reduce an optimal control problem to a two-point boundary value problem. Furthermore, for such functionals, the optimal controls can be obtained as explicit functions of the state and adjoint variables via the Pontryagin maximum principle. Although the financial cost of public health interventions to prevent epidemics most likely grows linearly with the magnitude of the corresponding control [24], we choose objective functionals with a quadratic dependence on the controls following the current trend of the literature [17, 18, 25–27].
4. The Pontryagin Maximum Principle
For control model (14) and the set of admissible controls , we can formulate two different optimal control problems using the objective functionals (16) and (17). Here, we prove the existence of solutions for these optimal control problems and we also obtain the optimality system in order to find a numerical approximation of the optimal controls.
Let us introduce the states and the controls and denote the right-hand side of the control model (14) as the vector function . For the objective functional (16), the optimal control problem of the HPV epidemic model is wheresubject to the dynamics of the control model (14) with given initial conditions.
We call a pair of states and controls satisfying both (14) and a feasible pair. heorem 4.1 in [28, Chapter III] ensures the existence of an optimal control and the corresponding solution for this problem. In particular, this existence theorem states that the following conditions are sufficient to guarantee the existence of an optimal control for (14):(H1) is continuous, and there exist positive constants and such that(a)(b) hold for all . Moreover, can be written as (H2)The set of admissible controls is closed and convex. Moreover, there is at least one feasible pair satisfying both (14) and (13)(H3) is convex on , and , , ,
Theorem 1. Consider the optimal control problem with control model (14) and cost functional (16). Then there exist such that .
Proof. Due to the a priori boundedness of the solutions of the control model (14) and the differentiability of the function , it follows that (a) and (b) in (H1) are ensured by suitable bounds on the partial derivatives of and on . Moreover, the state equations are linear with respect to the controls , and thus . The existence of a feasible pair is guaranteed by the Caratheodory theorem [29, pp. 182] for initial value problems. Moreover, for bounded controls on a finite time interval, is clearly closed and convex, and, hence, (H2) is satisfied as well. The integrand (19) of the cost functional is positive and quadratic with respect to the controls; therefore, is convex on . Furthermore, defining , we have , thus verifying (H3) and completing the proof.
Our next aim is to obtain the so-called optimality system in order to find numerically optimal controls , , , , . We define the HamiltonianBy the Pontryagin maximum principle [28], given the optimal control that minimizes the Hamiltonian (20) for almost all and the corresponding optimal solutions , there are adjoint variables that satisfy the systemwith the corresponding conditions for . From the optimality conditions, we obtain the following characterization for the optimal controls:
The proof of the existence of solutions for the optimal control problem with the objective functional (17) follows the same reasoning as in Theorem 1 and is, therefore, omitted. The technical details for the computation of the optimality system and characterization of the optimal control associated with objective functional (17) are provided in the appendix.
5. Results
We start by exploring the control model (14) numerically. The optimal control solutions can be obtained solving the optimality system. In order to obtain approximations of our optimality systems, we use the forward-backward sweep method (FBSM) from [30]. This iterative scheme begins with an initial guess on the controls; then, the state equations are solved forward in time. After that, since the transversality conditions for our problems are posed at the final time , the system for the adjoint variables is solved backward in time. The controls are updated using a convex combination of the previous controls and the new control obtained substituting the states and adjoints into its characterization. This process is repeated until convergence criteria are satisfied. For solving the state and adjoint systems, we use a fourth order Runge-Kutta scheme.
5.1. Model Parameters
Parameters values, units, and source of estimation are summarized in Table 1. We use these values unless otherwise is stated. Next, we outline their selection:(i)Vaccine efficacy, . Although clinical trials are still underway, there is evidence indicating efficacy of around for the nonavalent HPV vaccine Gardasil-9 [3]. We set .(ii)Vaccine protection, . The exact duration of vaccine protection is unknown, but clinical trials [19] have shown sustained efficacy for more than 5 years. We assume that the protection is from 5 to 15 years; thus years.(iii)Transmission rates, , , and . Transmission rates are the product of the average number of sexual contacts per unit of time and the probability of infection transfer per contact. We assume that both genders have between 1 and 5 sexual partners per year and per-act transmission probability range [3]. We set with . We also assume that the aware infected females take measures, such as condom use, that decrease the probability of infection transfer. Thus .(iv)Females’ infectious period, . The mean duration of a new genital HPV infection in women is 14.8 months for oncogenic types and 11.1 months for nononcogenic types [20]. Since the majority of the types targeted by Gardasil-9 are oncogenic, we assume that the female infectious period is between 9 months and 2 years; thus .(v)Male’s infectious period, . For men, most of the HPV infections clear in less than 12 months. Studies (see [21] and the references therein) report average clearance time of to months. We set year.(vi)Fraction of females that become aware of their infection, . Since is a fraction, .(vii)Periods of sexual activity, , and . We assume that the average duration of the sexual activity period for both, males and females, is between 15 and 50 years; therefore years.
Estimates of weight parameters associated with the objective functionals exhibit considerable variability in the literature. Thus, in the numerical simulations, several combinations of weight parameters are explored.
5.2. Uncontrolled System and Constant Control
In order to demonstrate the importance of the optimal controls, we firstly show the behavior of the model in the absence of control and under constant controls. The results are shown in Figure 1. It is easy to see that, without control, the system converges to an endemic equilibrium state (see Figure 1(a)). With constant controls (see Figure 1(b)) the infectious cases are reduced but the infection is not eradicated for this time interval. Besides, the number of vaccinated individuals tends to a very high proportion.
(a)
(b)
5.3. Optimal Vaccination and Screening
For the numerical computation of the optimal controls, we should complement the fixed parameters in Table 1 with values for the weights in the objective functionals. The cost of the nonavalent vaccine Gardasil-9 is approximately $134.26 per dose [18]. However, the associated weight parameters in functionals (16) and (17), namely, and , also include costs attributable to vaccine distribution, storage, and other administration-associated costs. It can be expected that administering the vaccine in school preadolescents is less expensive than vaccination of the sexually active adults; that is, [17]. Note also that the weights and are related to the cost of having infected individuals for this; they are different from weights , , which reflect the cost of the control interventions. Thus, the values of and and their rations to , , and are particularly important for the optimal controls. For example, if , that is, if treating infected individuals (including the cost of cancer induced by HPV) is assumed less costly than mass vaccination, then, in terms of minimizing cost, it can be logical to stop intervention programs. Therefore, in such a case, one can expect low vaccination and screening rates as the optimal policies. This consideration is confirmed by the results of the computations shown in Figure 2 where we plotted the optimal control corresponding to the cost functional (16) for
(a)
(b)
Please, observe that for this case the magnitude of the controls is very small (Figure 2(b)) and, consequently, the optimal controls are insufficient to reduce the spread of the infection (Figure 2(a)). Therefore, in what follows, we assume that the costs of having infected individuals are greater than the costs of administering the vaccine and screening programs. In particular, we assume that the values of weights and are two or three orders of magnitude bigger than the value of , , and . As a base case, we also set and . Before going any further to computations, we stress that, for problems of infectious disease control, the initial conditions can play an important role in the search of optimal controls. However, most studies overlook this issue. In this paper, we consider different initial conditions and explore how the optimal control varies depending on them. We choose three different sets of initial conditions to reflect precise infection levels, namely,(i)Low prevalence levels for both females and males (Figures 3(a), 3(b), 4(a), and 4(b))(ii)High prevalence for females and low prevalence for males (Figures 3(c), 3(d), 4(c), and 4(d))(iii)Low prevalence for women and high prevalence for males (Figures 3(e), 3(f), 4(e), and 4(f))
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
We start the numerical simulations considering the control model (14) with the objective functional (16) that aims to protect the vulnerable group. In other words, we seek to minimize the number of infections and cervical cancer cases in females. In this scenario, the optimal controls are characterized by system (22). In Figure 3, we show the numerical solution obtained via the FBSM. Figures 3(b), 3(d), and 3(f) illustrate the profiles of the optimal vaccination and screening rates, whereas Figures 3(a), 3(c), and 3(e) show the dynamics of the state variables. For the three initial conditions, we observe that and . This implies that the vaccination of females is more significant than increasing male’s immunization coverage. Nevertheless, we remark that, regardless of the independence of the objective functional (16) from male prevalence, the magnitude of the vaccination rates for males and females are nearly identical. Thus, to a certain extent, males' vaccination can be as important as females' vaccination. Furthermore, comparing Figures 3(b), 3(d), and 3(f), one can see that the optimal control profiles for the screening and vaccination rates are qualitatively and quantitatively similar. In other words, the optimal controls do not depend on the initial conditions. Please note that for all the scenarios explored, the optimal controls start with very high vaccination rates which quickly reduce the number of infected individuals to very low levels. These controls are then gradually reduced in value but are maintained enough time to achieve the ultimate eradication of the infection.
Our next task is to investigate the profiles of the optimal vaccination and screening rates with the cost functional (17) that aims to minimize the total prevalence at the final time . We compute the numerical solution of the two-point boundary value problem for the maximum principle via the FBSM (Figure 4). It is notable that also for the functional (17), the optimal vaccination rates for both sexes are remarkably similar for each of the initial conditions, respectively (see Figures 4(b), 4(d), and 4(f)). Even so, in contrast with the controls for functional (16), this time the controls do not start with high rates. This allows the number of infected individuals in the population to grow during the early phase of the epidemic. In spite of this, the optimal controls manage to eliminate the infection at the end of the time interval. Moreover, the infection levels in Figures 4(a), 4(c), and 4(e) are clearly higher than those in Figures 3(a), 3(c), and 3(e), respectively. Therefore, even though the controls for both functionals achieve the elimination of the infection, the controls for protection of the vulnerable group do it more efficiently. Nevertheless, there is considerably more vaccination in this case than for the case of eradication of the infection. Thus, the controls for functional (16) are presumably more expensive than the controls for functional (17).
Figure 5 presents a more direct comparison of relevant model outcomes under objective functionals (16) and (17). Figure 5(a) shows the total infected population as function of time under the optimal controls. Figure 5(b) shows the function which represents the cost of applying the optimal controls (without the cost of having infected individuals in the population). Figure 5(c) shows the reproduction number in the presence of vaccination and screening as function of time using the time series of the optimal controls. The blue dashed lines are the outcomes associated with functional (16) whose aim is to protect the vulnerable group. Meanwhile, the red solid lines are associated with functional (17) whose aim is to eradicate the infection.
(a)
(b)
(c)
From Figure 5 it is clear that the total prevalence is lower, for almost all times, when the objective is the protection of the vulnerable group (Figure 5(a)). Probably, this is a consequence of the huge difference in the magnitude of the vaccination rates for both functionals at the beginning of the time interval. Nevertheless, the cost of applying the control policies is less, on average, for the objective of eradication of the infection (Figure 5(b)). Yet, in the second half of the time interval, there is more vaccination occurring for the objective of eradication, resulting in greater costs for the controls associated with functional (17). This also explains why the value of the reproduction number is less for the objective of protection at the beginning of the time interval, but after a while, its value increases exceeding the value of the associated with the objective of eradication (Figure 5(c)). Finally, we stress that the cost of having infected individuals (which is of paramount importance) is not considered in Figure 5(b), therefore it is still uncertain which policy has the highest cost-effectiveness under more general hypotheses.
6. Discussion and Conclusion
In this paper, we investigate the optimal control policies against HPV infection. To address this problem, we introduced a two-sex compartmental epidemic model for the transmission dynamics of HPV in a heterosexual population. We incorporated five controls into this model which are able to act simultaneously. These five controls comprise all practical intervention policies used against HPV. Using the sex-specific infection transfers that we found, we computed the basic reproduction number in the presence of vaccination and screening. Using the van den Driessche theorem [23], we established the local stability of the disease-free equilibrium for values of less than one. This stability implies that, if initial infection levels are sufficiently low, maintaining ensures disease’s elimination and ultimate eradication of the virus. However, such a constant control scheme is usually a fragile strategy in a sense that it does not consider the changing dynamics of the disease’s spread and can be nonoptimal in terms of cost-efficacy.
A key issue in the implementation of HPV-vaccination is whether vaccinating both sexes is more effective than vaccinating the females only (as it is currently done in most countries). Moreover, although vaccinating individuals before their potential exposure to HPV infection can be a strategy with a positive outcome, vaccination of older individuals may be necessary if vaccine-induced protection is not lifelong. Here, we applied the optimal control theory to identify an effective distribution strategy for HPV-vaccines. We considered two different optimal control problems. Firstly, we studied a rather usual problem in the literature: minimizing the cumulative number of infected females and the cost of intervention policies over a given time interval. Such a problem corresponds to a policy that aimed at protection of the most vulnerable group of the population, that is, females at risk of developing invasive cancer. Secondly, we considered the problem of minimizing the total infection level in both sexes at the end of a given time interval at the lowest possible costs. This problem has a more ambitious objective, namely, the ultimate eradication of HPV. Such a goal is possible if the basic reproduction number is comparatively low. The possibility of total eradication of cervical cancer worldwide is currently discussed in the scientific community [31].
Results of our numerical simulations show that vaccination should be administered both prior and after the sexual debut and in both, females and males. For both functionals, the simulations support the intuitively expected outcome that females’ vaccination should be the priority; however, at the same time, we show that males’ vaccination is almost as beneficial and required. If the objective is to protect the vulnerable group, the numerical experiments imply that initially, higher vaccination rates should be applied which can gradually be reduced. This policy leads to a fast reduction in the number of infectious individuals. For the objective of eradication of the infection, the optimal controls begin at low rates allowing the prevalence to increase. Yet, the optimal vaccination and screening rates increase over time achieving the elimination of the infection at the end of the time interval.
We have to stress that our results differ from the majority of previous studies where the suggested preferred strategy was the vaccination of the preadolescent females prior to their sexual initiation. We believe that such a preference for the young females’ vaccination that prevails in the literature arose because the majority of authors did not consider the problem of HPV control as an optimal control problem. This also implies that the “young females only” vaccination is neither the optimal nor the most cost-effective policy. Furthermore, the efficacy of an HPV reduction program may also depend on the actual vaccine coverage. Our results were obtained for initially low vaccine coverage. This corresponds to a setting where HPV vaccination is just being introduced. Nevertheless, due to a possible nonlinear growth of the costs with respect to the fraction of the vaccinated individuals, even if the vaccination coverage of females is already high, it may still be cost-effective to include males into the vaccination program. Besides, there can be a group who is negligent or reluctant towards the vaccination; thus, a 100% vaccination coverage can be difficult to achieve even in one sex. Finally, we like to note that the impacts of vaccine coverage on the profiles of the optimal vaccination rates and the dependency of principal properties of the controls on the forms of objective functionals remain open questions and deserve further study.
Appendix
Here, we compute the adjoint system and the characterization of the optimal controls considering the objective functional (17): In this case, the optimal control problem of the HPV epidemic model issubject to the dynamics of the control model (14) with given initial conditions.
The Pontryagin maximum principle converts the optimal control problem (A.2) into a problem of minimizing pointwise the Hamiltonianwith respect to the controls. Here, for are adjoint variables that satisfy the system of differential equationswith transversality conditions for , and for . For this control problem, applying the maximum principle from [32], we obtain the following characterization for the optimal controls:
Data Availability
The software and data of computation used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Fernando Saldaña acknowledges the support by CONACyT, Mexico, under grant number 331294. Andrei Korobeinikov is supported by Spain’s Ministerio de Ciencia, Innovación y Universidades grant MTM2015-71509-C2-1-R.