For each analysis presented below, we assume that a couple has generated, by IVF, $n$ viable embryos such that each embryo, if implanted, would have led to a live birth. We focus on a single complex disease, and assume that the corresponding PRS has been computed for each embryo. Given the PRSs of the $n$ embryos, a single embryo is selected for implantation based on a selection strategy.

The first strategy we consider is aimed only at avoiding high-risk embryos, consistent with studies of the potential clinical utility of PRSs in adults (Chatterjee et al., 2016; Dai et al., 2019; Gibson, 2019; Khera et al., 2018; Mars et al., 2020; Mavaddat et al., 2019; Torkamani et al., 2018). For example, the first case report presented on PES described the identification and exclusion of embryos with extremely high (top 2-percentiles) PRS (Treff et al., 2019a). We term this strategy ‘high-risk exclusion’ (HRE: Figure 1A, upper panel). Under HRE, after high-risk embryos are set aside, an embryo is randomly selected for implantation among the remaining available embryos. (In the case that all embryos are high-risk, we assume a random embryo is selected among them.)

Figure 1

Download asset Open asset

An alternative selection strategy is to use the embryo with the lowest PRS. Ranking and prioritizing embryos for implantation based on morphology is common in current IVF practice (Bormann et al., 2020; Montag et al., 2013; Rhenman et al., 2015). If ranking is instead based on a disease PRS, the embryo with the lowest PRS could be selected, without any recourse to high-risk PRS thresholds. Such an approach was suggested by another recent publication from the same company (based on a multi-disease index), but outcomes were only examined in the context of sibling pairs (Treff et al., 2020). We term the implantation of the embryo with the lowest PRS as ‘lowest-risk-prioritization’ (LRP; Figure 1A, lower panel).

In the following, we describe the theoretical risk reduction that can be achieved under these selection strategies. Our statistical approach is based on the liability threshold model (LTM; Falconer, 1967). The LTM represents disease risk as a continuous liability, comprising genetic and environmental risk factors, under the assumption that individuals with liability exceeding a threshold are affected. The liability threshold model has been shown to be consistent with data from family-based transmission studies (Wray and Goddard, 2010) and GWAS data (Visscher and Wray, 2015). Consequently, we define the disease risk of a given embryo probabilistically, as the chance that, given its PRS, its liability will cross the threshold at any point after birth (Figure 1B).

We use the following notation. We define the predictive power of a PRS as the proportion of variance in the liability of the disease explained by the score (Dudbridge, 2013), and denote it as ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2}$. We quantify the outcome of PES in two ways: the relative risk reduction (RRR) is defined as $\text{RRR}=\frac{K-P\left(\text{disease}\right)}{K}=1-\frac{P\left(\text{disease}\right)}{K}$, where $K$ is the disease prevalence and $P\left(disease\right)$ is the probability of the selected embryo to be affected; the absolute risk reduction (ARR) is defined as $K-P\left(disease\right)$. For example, if a disease has prevalence of 5% and the selected embryo has a probability of 3% to be affected, the RRR is 40%, and the ARR is 2% points. We computed the RRR and ARR analytically under each selection strategy, and for various values for the disease prevalence, the strength of the PRS, embryo exclusion thresholds, and other parameters. The mathematical basis of the calculations is summarized in Materials and methods, and detailed in the Appendix.

In Figure 2 (upper row), we show the relative risk reduction achievable under the HRE strategy with $n=5$ embryos. Under the 2-percentile threshold (straight black lines), the reduction in risk is limited: the RRR is <10% in all scenarios where ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2\le 0.1}$. Currently, ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2\approx 0.1}$ (on the liability scale) is the upper limit of the predictive power of PRSs for most complex diseases (Lambert et al., 2021), with the exception of a few disorders with large-effect common variants (such as Alzheimer’s disease or type 1 diabetes) (Sharp et al., 2019; Zhang et al., 2020). In the future, more accurate PRSs are expected. However, the common-variant SNP heritability is at most ≈30% even for the most heritable diseases such as schizophrenia and celiac disease (Holland et al., 2020; Zhang et al., 2018), and it was recently suggested that ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2=0.3}$ is the maximal realistic value for the foreseeable future (Wray et al., 2021). At this value, relative risk reduction would be 20% for $K=0.01$, 9% for $K=0.05$, and 3% for $K=0.2$. These gains achieved with HRE are small because the overwhelming majority of affected individuals do not have extreme scores (Murray et al., 2021; Wald and Old, 2019).

Figure 2 with 3 supplements see all

Download asset Open asset

Risk reduction increases as the threshold for exclusion is expanded to include the top quartile of scores, and then reaches a maximum at ≈25-50% under a range of prevalence and ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2}$ values. For all these simulations, we set the number of available (testable) embryos to $n=5$ (Dahdouh, 2021; Sunkara et al., 2011), although we acknowledge that the number of viable embryos may be much lower for many couples seeking IVF services for infertility (Smith et al., 2015). Simulations show that these estimates do not change much with increasing the number of embryos (see Figure 2—figure supplement 1). This holds especially at more extreme threshold values, since most batches of $n$ embryos will not contain any embryos with a PRS within, for example, the top 2-percentiles.

It should be noted that the relative risk reduction does not increase monotonically under HRE. Under our definition, whenever all embryos are high risk, an embryo is selected at random. Thus, at the extreme case when all embryos (i.e. top 100%) are designated as high risk, an embryo is selected at random at all times, and the relative risk reduction reduces to zero. We chose this definition of the HRE strategy because it does not involve ranking of the embryos. However, we can also consider an alternative strategy: if all embryos are high risk, the embryo with the lowest PRS is selected. Here, the RRR is expected to increase when increasing the threshold and designating more embryos as high risk, which we confirm in Figure 2—figure supplement 2. When the threshold is at 100% (all embryos are high risk), this alternative strategy (which we do not further consider) reduces to the lowest-risk prioritization strategy, which we study next.

The HRE strategy treats all non-high-risk embryos equally. In practice, we expect most, or even all, embryos to be designated as non-high-risk, given the recent focus on the top PRS percentiles in the literature (e.g. Khera et al., 2018). However, as we have seen, this strategy leads to very little risk reduction. In Figure 2 (lower panels), we show the expected RRR for the lowest-risk prioritization strategy, under which we prioritize for implantation the embryo with the lowest PRS, regardless of any PRS cutoff. Indeed, under the LRP strategy, risk reductions are substantially greater than in HRE. For example, with $n=5$ available embryos, RRR>20% across the entire range of prevalence and ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2}$ parameters considered, and can reach ≈50% for $K\le 5\%$ and ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2=0.1}$, and even ≈80% for $K=1\%$ and ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2=0.3}$. While RRR continues to increase as the number of available embryos increases, the gains are quickly diminishing after $n=5$. On the other hand, Figure 2 also demonstrates that RRR drops steeply if the number of embryos falls below $n=5$, although the lower bound for RRR when just two embryos are available (≈20% for many scenarios) is still comparable to the upper bound of the HRE strategy for a greater number of embryos.

Our results demonstrate that, contrary to our previous study reporting only small effects of PES for quantitative traits (Karavani et al., 2019), PES can generate substantial relative risk reductions for diseases under the LRP strategy. To understand the relation between continuous and binary traits, consider an example involving IQ. Our estimate for the mean gain in IQ that could be achieved by selecting the embryo with the highest IQ polygenic score is approximately ≈2.5 IQ points (Karavani et al., 2019). Now assume that individuals with IQ<70 (2 SDs below the mean) are considered 'affected' according to a dichotomized trait of 'cognitive impairment'. Among individuals with IQ<70, the proportion of individuals with IQ in the range [67.5,70] is 33.5% (assuming a normal distribution). A gain of 2.5 points would shift such offspring beyond the threshold for 'cognitive impairment', resulting in a corresponding 33.5% reduction in risk of being 'affected'. (Note that the above explanation is intended to provide an intuition and ignores any variability in the gain.) Figure 2—figure supplement 3 utilizes statistical modeling (with ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2}$ derived from recent GWAS for intelligence [Savage et al., 2018]) to demonstrate that substantial risk reductions can be achieved for a dichotomized trait, including when selecting out of just three embryos (panel (A)). Panel (B) extends these results to data for LDL cholesterol (with ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2}$ derived from Weissbrod et al., 2021); given $n=5$ embryos and the currently available PRS for LDL-C levels, risk reductions for 'high cholesterol' range from 40 to 60%, depending on the LDL level used to define the categorical trait. Thus, while implanting the embryo with the most favorable PRS is expected to result in very modest gains in an underlying quantitative trait, it is at the same time effective in avoiding embryos at the unfavorable tail of the trait.

We next examined the effects of parental PRSs on the achievable risk reduction (Materials and methods, see also the Appendix), given that families with high genetic risk for a given disease may be more likely to seek PES. Figure 3 demonstrates that, as expected, the HRE strategy shows greater relative risk reduction as parental PRS increases, in particular when excluding only very high-scoring embryos. This result follows directly from the fact that, on average, offspring will tend to have PRS scores near the mid-parental PRS value. In contrast, the relative RR (although not the absolute RR; see next section) for the LRP strategy somewhat declines as parental PRSs increase. Nevertheless, the RRR for the LRP strategy remains greater than that for the HRE strategy across all parameters (as expected by the definitions of these strategies).

Figure 3 with 2 supplements see all

Download asset Open asset

It is also conceivable that families may be more likely to seek PES when one or both prospective parents is affected by a given disease. In Figure 3—figure supplement 1, we plot the RRR under the HRE and LRP strategies given that the parents are both healthy, both affected, or one of each (where we fixed the prevalence $K=5\%$ and the heritability to $h^2=40\%$). The figure illustrates that parental disease status has relatively little impact on the expected RRR (especially in comparison to the changes under HRE when conditioning on the actual parental PRSs). This is because, as long as ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2\ll 1,}$ parental disease does not necessarily provide much information about parental PRS, and thus does not strongly constrain the number of risk alleles available to each embryo.

The above results were presented in terms of relative risk reductions. However, Figure 3—figure supplement 1 also shows the baseline risk of an embryo of parents with a given disease status. For example, when one of the parents is affected, selecting the lowest risk embryo out of $n=5$ (for a realistic ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2=0.1}$) reduces the risk from 10.0% to only 5.8%, thus nearly restoring the risk of the future child to the population prevalence (5%). More generally, we plot the absolute risk reduction (ARR) under the HRE and LRP strategies in Figure 3—figure supplement 2 for a few values of the parental PRSs. Notably, while RRRs under the LRP strategy somewhat decrease with increasing parental PRSs, the ARRs substantially increase, in accordance with an expectation that PES in higher risk parents should eliminate more disease cases.

The clinical interpretation of these absolute risk changes will vary based on the population prevalence of the disorder (or the baseline risk of specific parents), and can offer a very different perspective on the magnitude of the effects (Gordis, 2014; Lázaro-Muñoz et al., 2021; Murray et al., 2021). In particular, for a rare disease, large relative risk reductions may result in very small changes in absolute risk. As an example, schizophrenia is a highly heritable (Sullivan et al., 2003) serious mental illness with prevalence of at most 1% (Perälä et al., 2007). The most recent large-scale GWAS meta-analysis for schizophrenia (Ripke et al., 2020) has reported that a PRS accounts for approximately 8% of the variance on the liability scale. Our model shows that a 52% RRR is attainable using the LRP strategy with $n=5$ embryos. However, this translates to only ≈0.5 percentage points reduction on the absolute scale: a randomly-selected embryo would have a 99% chance of not developing schizophrenia, compared to a 99.5% chance for an embryo selected according to LRP. In the case of a more common disease such as type 2 diabetes, with a lifetime prevalence in excess of 20% in the United States (Geiss et al., 2014), the RRR with $n=5$ embryos (if the full SNP heritability of 17% [Zhang et al., 2018] were achieved) is 43%, which would correspond to >8 percentage points reduction in absolute risk.

The results depicted in Figure 2 describe the average risk reduction across the population, whereas the results in Figure 3 demonstrate results for specific combinations of parental risk scores. However, it remains unclear whether the large average risk reductions observed under the LRP strategy are driven by only a small proportion of couples. More generally, we would like to fully characterize the dependence of the risk reduction on parental PRSs, which could be of interest to physicians and couples in real-world settings.

To address these questions, we define a new risk reduction index, which we term the per-couple relative risk reduction, or pcRRR. Informally, the pcRRR is the relative risk reduction conditional of the PRSs of the couple. Mathematically, $\text{pcRRR(couple)}=1-\frac{P_s\left(\text{disease}|\text{couple}\right)}{P_r\left(\text{disease}|\text{couple}\right)}$. Here, ${\displaystyle P_s\left(\text{disease}|\text{couple}\right)}$ is the probability that the (PRS-based) selected embryo is affected given the PRSs of the couple, and ${\displaystyle P_r\left(\text{disease}|\text{couple}\right)}$ is similarly defined for a randomly selected embryo. Conveniently, the pcRRR depends only on the average of the maternal and paternal PRSs, which we denote as $c$. We calculated $\text{pcRRR}\left(c\right)$ analytically under the LRP strategy (the Appendix), as well as computed the distribution of $\text{pcRRR}\left(c\right)$ across all couples in the population.

We show the distribution of $\text{pcRRR}\left(c\right)$ in Figure 4, panels (A)-(C). The results demonstrate that the pcRRR is relatively narrowly distributed around its mean, for all values of the prevalence ($K$) considered. The distribution becomes somewhat wider (and left-tailed) for the most extreme ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2}$ (0.3). Thus, the population-averaged RRRs are not driven by a small proportion of the couples. In agreement, the pcRRR depends only weakly on the average parental PRS, as can be seen in panels (D)-(F).

Figure 4

Download asset Open asset

We note that the per-couple relative risk reduction is also an average, over all possible batches of $n$ embryos of the couple. One may thus ask what is the distribution of possible RRRs across these batches. We provide a short discussion in the Appendix (Section 5.3).

Polygenic risk scores are often correlated across diseases (Watanabe et al., 2019; Zheng et al., 2017). Therefore, selecting based on the PRS of one disease may increase or decrease risk for other diseases. While a full analysis of screening for multiple diseases is left for future work, our simulation framework allows us to investigate the potential harmful effects of prioritizing embryos for one disease, in case that disease is negatively correlated with another disease (see the Appendix). We considered genetic correlations between diseases taking the values $\rho =(-0.05,-0.1,-0.15,-0.2,-0.3)$. [The most negative correlation between two diseases reported in LDHub (https://ldsc.broadinstitute.org/ldhub/) is −0.3, occurring between ulcerative colitis and chronic kidney disease (Zheng et al., 2017).] In general, negative correlations between diseases are uncommon, and when they occur, typical correlations are about −0.1.

Figure 5 shows the simulated risk reduction for the target disease and the risk increase for the correlated disease, across different values of $\rho$ and for three values of the prevalence $K$ (panels (A)-(C); assumed equal for the two diseases), all under the LRP strategy. In all panels, we used ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2=0.1}$ for both diseases. The relative risk reduction for the target disease is, as expected, always higher in absolute value than the risk increase of the correlated disease. For typical values of $\rho =-0.1$ and $n=5$, the relative increase in risk of the correlated disease is relatively small, at ≈6% for $K\le 0.05$ and ≈3.5% for $K=0.2$. However, for strong negative correlation ($\rho =-0.3$) the increase in risk can reach 22%, 16%, or 11% for ${\displaystyle K=0.01,0.05}$ and ${\displaystyle 0.2}$, respectively. Thus, care must be taken in the unique setting when the target disease is strongly negatively correlated with another disease.

Figure 5

Download asset Open asset

Our analysis so far has been limited to mathematical analysis and simulations based on a statistical model. In principle, it would be desirable to compare our predictions to results based on real data. However, clearly, no real genomic and phenotypic data exist that would correspond to our setting, nor could such data be ethically or practically generated. Thus, we resort to a ‘hybrid’ approach, in which we simulate the genomes of embryos based on real genomic data from case-control studies. This approach is similar to the one we have previously used for studying polygenic embryo screening for traits (Karavani et al., 2019).

Briefly, our approach is as follows. We consider separately two diseases with somewhat differing genetic architecture: schizophrenia, which is amongst the most polygenic complex diseases, with no common loci of high effect size, and Crohn’s disease, which is estimated to be less polygenic, and has several common loci with much larger effects than those found in schizophrenia (O'Connor et al., 2019). For each disease, we used genomes of unrelated individuals drawn from case-control studies. For schizophrenia, we used ≈900 cases and ≈1600 controls of Ashkenazi Jewish ancestry, while for Crohn’s, we used ≈150 cases and ≈100 controls of European ancestry. We then generated 'virtual couples' by randomly mating pairs of individuals, regardless of sex. For each couple, we simulate the genomes of $n$ hypothetical embryos, based on the laws of Mendelian inheritance and by randomly placing crossovers according to genetic map distances. In parallel, we used the 'parental' genomes to learn a logistic regression model that predicts the disease risk given a PRS computed based on existing summary statistics. We then computed the PRS of each simulated embryo, and predicted the risk that embryo to be affected. Finally, we compared the risk of disease between a population in which one embryo per couple is selected at random, vs. a population in which one embryo is selected based on its PRS. For complete details, see Materials and methods.

In Figure 6, we plot the results for the relative risk reduction for schizophrenia (panels (A) and (B)) and Crohn’s disease (panels (C) and (D)). For each disease, we consider both the HRE and LRP strategies. The analytical predictions closely match the empirical risk reductions generated in the simulations, except for a slight overestimation of the RRR under the LRP strategy. Nevertheless, for both schizophrenia and Crohn’s disease, we empirically observe that RRRs as high as ≈45% are achievable with $n=5$ embryos. In contrast, under the HRE strategy and when excluding embryos at the top 2% risk percentiles, risk reductions are very small, in agreement with the theoretical predictions. These results thus provide support to the robustness of our statistical model.

Figure 6 with 1 supplement see all

Download asset Open asset

To further investigate the assumptions of our model, we test in Figure 6—figure supplement 1 two intermediate predictions. The first is that the variance of the PRSs of embryos of a given couple should not depend on the average parental PRS. This is indeed the case (panels (A) and (C)), with the only exception of an uptick of the variance at very low parental PRSs for schizophrenia. The second prediction is that the variance across embryos is half of the variance in the parental population. The empirical results again show reasonable agreement with the theoretical prediction (panels (B) and (D)). The empirical variance (averaged across couples) was slightly lower than expected (by ≈4% for schizophrenia and ≈14% for Crohn’s), which may explain our slight overestimation of the expected RRR under the LRP strategy.

In this section, we provide a brief overview of our model and derivations, with complete details appearing in the Appendix.

Our model is follows. We write the polygenic risk scores of a batch of $n$ IVF embryos as $\left(s_1,\dots ,s_n\right)$, and generate the scores as $s_i=x_i+c$. The $\left(x_1,\dots ,x_n\right)$ are embryo-specific independent random variables with distribution ${\displaystyle x_i\sim N\left(0,r_{\mathrm{p}\mathrm{s}}^2/2\right)}$, ${\displaystyle r_{\mathrm{p}\mathrm{s}}^2}$ is the proportion of variance in liability explained by the score, and $c$ is a shared component with distribution ${\displaystyle c\sim N\left(0,r_{\mathrm{p}\mathrm{s}}^2/2\right)}$, also representing the average of the maternal and paternal scores.

In each batch, an embryo is selected according to the selection strategy. Under high-risk exclusion, we select a random embryo with score ${\displaystyle s<z_q{r}_{\mathrm{p}\mathrm{s}}}$, where $z_q$ is the $\left(1-q\right)$-quantile of the standard normal distribution. If no such embryo exists, we select a random embryo, but we also studied the strategy when in such a case, the lowest scoring embryo is selected. Under lowest-risk prioritization, we select the embryo with the lowest value of $s$. We computed the liability of the selected embryo as $y=s+e$, where ${\displaystyle e\sim N\left(0,1-r_{\mathrm{p}\mathrm{s}}^2\right)}$. We designate the embryo as affected if ${\displaystyle y>z_K}$, where $z_K$ is the $(1-K)$ -quantile of the standard normal distribution and $K$ is the disease prevalence. In the simulations, we computed the disease probability (for each parameter setting) as the fraction of batches (out of 10⁶ repeats) in which the selected embryo was affected. We also simulated the score and disease status of a second disease, which is not used for selecting the embryo, but may be negatively correlated with the target disease.

We computed the disease probability analytically using the following approaches. We first computed the distribution of the score of the selected embryo. For lowest-risk prioritization, we used the theory of order statistics. For high-risk exclusion, we first conditioned on the shared component $c$, and then studied separately the case when all embryos are high-risk (i.e. have score ${\displaystyle s>z_q{r}_{\mathrm{p}\mathrm{s}}}$), in which the distribution of the unique component of the selected embryo ($x$) is a normal variable truncated from below at ${\displaystyle z_q{r}_{\mathrm{p}\mathrm{s}}-c}$, and the case when at least one embryo has score ${\displaystyle s<z_q{r}_{\mathrm{p}\mathrm{s}}}$, in which $x$ is a normal variable truncated from above. We then integrated over the non-score liability components (and over $c$ in some of the settings) in order to obtain the probability of being affected. We solved the integrals in the final expressions numerically in R.

We computed the risk reduction based on the ratio between the risk of a child of a random couple when the embryo was selected by PRS and the population prevalence. We also provide explicit results for the case when the average parental PRS $c$ is known. These expressions allowed us to compute the distribution of risk reductions per-couple. Finally, when conditioning on the parental disease status, we integrated the disease probability of the selected embryo over the posterior distribution of the parental score and non-score genetic components. For full details and for an additional discussion of previous work and limitations, see the Appendix. R code is available at: https://github.com/scarmi/embryo_selection (copy archived at swh:1:rev:4cdc572582deb9b745e6844d96e0344914f4595e, Carmi, 2021) and https://github.com/dbackenroth/embryo_selection (copy archived at swh:1:rev:c65bf082fcb28434c271260560c4a4450dad76a3,; Backenroth, 2021).

Our main analysis has been limited to mathematical modeling of polygenic scores and their relation to disease risk. For obvious ethical and practical reasons, we could not validate our modeling predictions with actual experiments. Nevertheless, we could perform realistic simulations based on genomes from case-control studies, similarly to our previous work (Karavani et al., 2019). Our approach is as follows. We consider, separately, two diseases: schizophrenia and Crohn’s. For schizophrenia, we use ≈900 cases and ≈1600 controls of Ashkenazi Jewish ancestry, while for Crohn’s, we use ≈150 cases and ≈100 controls from the New York area. For each disease, we use these individuals, who are unrelated, to generate 'virtual couples' by randomly mating pairs of individuals. For each such 'couple', we simulate the genomes of $n$ hypothetical embryos, based on the laws of Mendelian inheritance and by randomly placing crossovers according to genetic map distances. In parallel, we use the same genomes to derive a logistic regression model that predicts the risk of disease given a PRS computed from the most recently available summary statistics (based on datasets not including the samples in our test cohorts). We then compute the PRS of each simulated embryo, and predict the risk of disease of that embryo. We finally compare the risk of disease between one randomly selected embryo per couple vs one embryo selected based on PRS. In the paragraphs below, we provide additional details.

The liability threshold model (LTM) is a classic model in quantitative genetics (Dempster and Lerner, 1950; Falconer, 1965; Lynch and Walsh, 1998) and is also commonly used to analyze modern data (e.g. Wray and Goddard, 2010; So et al., 2011; Lee et al., 2011; Lee et al., 2012; Do et al., 2012; Hayeck et al., 2017; Weissbrod et al., 2018; Hujoel et al., 2020). Under the LTM, a disease has an underlying ‘liability’, which is normally distributed in the population, and is the sum of two components: genetic and non-genetic (the environment). Further, the LTM assumes an ‘infinitesimal’, or ‘polygenic’ genetic basis, under which a very large number of genetic variants of small effect combine to form the genetic component. An individual is affected if his/her total liability (genetic + environmental) exceeds a threshold.

Mathematically, if we denote the liability as $y$, the LTM can be written as

where $y\sim N(0,1)$ is a standard normal variable, $g\sim N(0,h^2)$ is the genetic component, with variance equals to the heritability $h^2$, and $ϵ\sim N(0,1-h^2)$ is the non-genetic component. In practice, we cannot measure the genetic component, but only estimate it imprecisely with a polygenic risk score, denoted $s$. Following previous work (So et al., 2011; Do et al., 2012; Lee et al., 2012; Treff et al., 2019a; Karavani et al., 2019), we assume that the LTM can be written, similarly to Equation (1), as

where $y\sim N(0,1)$ as above, $s\sim N(0,r_{\text{ps}}^2)$, where $r_{\text{ps}}^2$ is the proportion of the variance in liability explained by the score, and $e\sim N(0,1-r_{\text{ps}}^2)$ is the residual of the regression of the liability on $s$ (and is uncorrelated with $s$), representing environmental effects as well as genetic factors not accounted for by the score.

An individual is affected whenever his/her liability exceeds a threshold. The threshold is selected such that the proportion of affected individuals is equal to the prevalence $K$, that is, it is equal to $z_K$, the $(1-K)$-quantile of a standard normal variable. Thus,

The model is illustrated in Figure 1B of the main text.

Consider the polygenic risk scores (for a disease of interest) of $n$ IVF embryos of given parents. We assume no information is known about the parents, or, in other words, that the parents are randomly and independently drawn from the population. The scores of the embryos have a multivariate normal distribution,

 where the means form a vector ${\mathrm{}}_n$ of $n$ zeros, and the $n\times n$ covariance matrix is

The diagonal elements of the matrix are simply the variances of the individual scores of each embryo. The off-diagonal elements represent the covariance between the scores of the embryos, who are genetically siblings. Based on standard quantitative genetic theory (Lynch and Walsh, 1998) (see also our previous paper [Karavani et al., 2019]), the covariance between the scores of two siblings is $\text{Cov}(s_i,s_j)=\frac{1}{2}\text{Var}\left(s\right)$, and hence the off-diagonal elements follow. [The non-score components (the $e$ terms in Equation (2)) are also correlated. The correlation between the genetic components of $e$ is modeled in Section 6. Modeling the correlation between the environmental components was unnecessary in this paper – see Section 10].

As we showed in our previous work (Karavani et al., 2019), the scores can be written as a sum of two independent multivariate normal variables, $\mathit{𝒔}=\mathit{𝒙}+\mathit{𝒄}$, with

where ${\mathrm{}}_n$ is a vector of zeros of length $n$, ${\mathit{𝑰}}_n$ is the $n\times n$ identity matrix, and ${\mathit{𝑱}}_n$ is the $n\times n$ matrix of all ones. The x_i’s and c_i’s have the same marginal distribution, namely normal with zero mean and variance $r_{\text{ps}}^2/2$ each. However, the x_i’s are independent, whereas $\mathit{𝒄}$ has a constant covariance matrix, which means that the c_i’s are $n$ identical copies of the same random variable,

Thus, for each embryo

$i=1,\mathrm{\dots },n$, 

We assume next that we select for implantation the embryo with the lowest polygenic risk score for the disease of interest. Our goal will be to calculate the probability of that embryo to be affected. Since $s_i=x_i+c$, the score of the selected embryo satisfies 

where we defined $x_{\min }=\min (x_1,\mathrm{\dots },x_n)$. Denote by $i^{*}$ the index of the selected embryo ($x_{i^{*}}=x_{\min }$). The liability of the embryo with the lowest risk is thus

where e_i is the non-score component of embryo $i$, and $\tilde{e}=c+e_{i^{*}}$. We have,

Therefore, the liability of the selected embryo can be written as a sum of two (independent) variables: x_min, which is the minimum of $n$ independent (zero mean) normal variables with variance $r_{\text{ps}}^2/2$ each; and $\tilde{e}$, which is a normal variable with (zero mean and) variance $1-r_{\text{ps}}^2/2$.

The distribution of x_min can be computed based on the theory of order statistics,

In the above equation, the minimum of $n$ variables is greater than $t$ if and only if all variables are greater than $t$. The distribution of each $x$ is normal with zero mean and variance $r_{\text{ps}}^2/2$, and hence $P(x>t)=1-\mathrm{\Phi }\left(\frac{t}{r_{\text{ps}}/\sqrt{2}}\right)$, where $\mathrm{\Phi }(\cdot )$ is the cumulative probability distribution (CDF) of a standard normal variable.

We can now compute the probability of the selected embryo to be affected by demanding that the total liability is greater than the threshold $z_K$. Denote the probability of disease as $P_s(\text{disease})$ ($s$ stands for selected). Conditional on $\tilde{e}$,

where in the fourth line, we used Equation (18). Next, denote by $f(\tilde{e})$ the density of $\tilde{e}$, and by $\varphi (\cdot )$ the probability density function of a standard normal variable. Given that $\tilde{e}\sim N(0,1-r_{\text{ps}}^2/2)$,

In the third line, we changed variables: $t=\tilde{e}/\sqrt{1-r_{\text{ps}}^2/2}$. Equation (20) is our final expression for the probability of the embryo with the lowest score to be affected.

We now consider the selection strategy in which the implanted embryo is selected at random, as long as its risk score is not particularly high. Specifically, we assume that whenever possible, embryos at the top $q$ risk percentiles are excluded. When all embryos have high risk, we assume that a random embryo is selected. Let z_q be the $(1-q)$-quantile of the standard normal distribution. The variance of the score is $r_{\text{ps}}^2$, and therefore, the score of the selected embryo must be lower than $z_q{r}_{\text{ps}}$.

To compute the disease risk in this case, we first condition on the shared, family-specific component $c$. We later integrate over $c$ to derive the risk across the population. Denote by x_s the value of $x$ for the selected embryo, and for the moment, also condition on x_s. We have,

To obtain $P_s(\text{disease}|c)$, we need to integrate over $f(x_s)$, the density of x_s. In fact, $f(x_s)$ is a mixture of two distributions, depending on whether or not all embryos were high risk. Denote by $H$ the event that all embryos are high risk, and let us first compute the probability of $H$. Recall that given $c$, the scores of all embryos, $s_i=x_i+c$, are independent. The event $H$ is equivalent to the intersection of the independent events $\{s_i>z_q{r}_{\text{ps}}\}$ for $i=1,\mathrm{\dots },n$. Thus, recalling that $x_i\sim N(0,r_{\text{ps}}^2/2)$,

Given $H$, we know that all scores were higher than the cutoff, i.e., that $x_i>z_q{r}_{\text{ps}}-c$ for all $i=1,\mathrm{\dots },n$. An embryo is then selected at random. Thus, x_s, the value of $x$ of the selected embryo, is a realization of a normal random variable truncated from below. Specifically, if $f_x(\cdot )$ is the unconditional density of $x$, then for $x_s>z_q{r}_{\text{ps}}-c$,

In the case $H$ did not occur, we select an embryo at random among embryos with score $s_i<z_q{r}_{\text{ps}}$, that is, $x_i<z_q{r}_{\text{ps}}-c$. The density of x_s is again, analogously to the above case, a realization of a normal random variable, but this time truncated from above. For

$x_s<z_q{r}_{\text{ps}}-c$,

Using these results, we can write the density of x_s when conditioning only on $c$,

We can now integrate over all x_s, still conditioning on $c$, and using Equation (24) and some algebra,

where we defined

Equation (29) provides an expression for the probability of a disease given the mean parental score $c$.

Finally, we can integrate over all $c$ in order to obtain the probability of disease in the population. Recalling that $c\sim N(0,r_{\text{ps}}^2/2)$ and denoting its density as $f(c)$, and again after some algebra,

where we defined

and $\eta (t,\cdot )$ was defined in Equation (30) above. Equation (31) is our final expression for the probability of an embryo to be affected after being selected randomly among non-high-risk embryos.

We define the relative risk reduction (RRR) as follows. We are given the prevalence $K$ and the probability of the selected embryo to be affected $P_s(\text{disease})$ (averaged over the population). Then,

The absolute risk reduction (ARR) is similarly defined as $K-P_s(\text{disease})$. For example, if a disease has prevalence of 5% and an embryo selected based on PRS has an average probability of 3% to be affected, the relative risk reduction is 40%, while the absolute risk reduction is 2% points.

To use Equation (33), $P_s(\text{disease})$ is given by Equation (20) for the lowest-risk prioritization strategy, and by Equation (31) for the high-risk exclusion strategy. We solve the integrals in these equations numerically in R using the function integrate (see Section 11).

In the following, we compute the relative risk reduction when the disease status of the parents is given.