A logical network-based drug-screening platform for Alzheimer’s disease representing pathological features of human brain organoids

Developing effective drugs for Alzheimer’s disease (AD), the most common cause of dementia, has been difficult because of complicated pathogenesis. Here, we report an efficient, network-based drug-screening platform developed by integrating mathematical modeling and the pathological features of AD with human iPSC-derived cerebral organoids (iCOs), including CRISPR-Cas9-edited isogenic lines. We use 1300 organoids from 11 participants to build a high-content screening (HCS) system and test blood–brain barrier-permeable FDA-approved drugs. Our study provides a strategy for precision medicine through the convergence of mathematical modeling and a miniature pathological brain model using iCOs.

* Research reproducibility of the modeling results: While the Boolean functions of the model are available as plain text in the Supplementary Excel file, the model should be shared in the SBMLqual format (PMID: 24321545) to make it easier for others to reproduce the presented research, and to further expand the model. At a minimum, this SBML file should be available as part of this manuscript. Better would be making the model available in any or all of existing Boolean model repositories (e.g., BioModels (which has been supporting Boolean models for a while now; PMID: 29106614 ), GINsim.org, Cell Collective; PMID: 22871178). * Related to reproducibility, I don't believe the methods for the modeling component are comprehensive enough. For example, what initial state(s) did the authors use for the quality control input-output analyses ( Fig. 5 and Fig. 5-1), and the actual predictions? Did the authors perform initial condition sensitivity analysis? * For the oxidative stress input-output analyses: what activity levels were used for other external variables (e.g., Reelin, APOE4, etc.)? Did the authors do these analyses for multiple external conditions?
By the way, the idea of phenotype score is interesting to estimate the pathological level, but there are several concerns. Authors choose key proteins and phenotypes such as Aβ, p-tau, synapse loss, apoptosis, and autophagy for calculation of phenotype score, but how to choose these key proteins and phenotypes? Phenotype score looks work well but p-tau proteins have double impact on the score because p-tau activates synapse loss. They developed a network model but they did not consider these kinds of network effects in the calculation of phenotype score. If phenotype score well works to estimate the pathological level, why did they need to develop a network model and conduct a perturbation analysis? Why don't authors directly conduct the phenotype score calculation? Basic idea of this manuscript is a novel and interesting in developing drug-screening for Alzheimer's disease, but they need to clearly describe their methods and results.
Reviewer #3: Remarks to the Author: This is an interesting paper in which Park et al. first generate a cerebral organoid (CO) model of Alzheimer's disease (AD) through two complementary appraoches: (1) using patient derived iPSCs (PiB+ iCOs) and (2) using isogenic iPSC lines with a CRISPR-edited apoE ε4 allele (E4iso iCOs).
They do an in-depth characterization of their AD organoids to show that they recapitulate some key pathological phenotypes of the disease, which are absent in other in-vitro models. Finally, they use computational modeling of known AD pathways to narrow down a list of candidate targets, which are selected using mathematical modelling, for testing FDA-approved drugs. They attempt to validate these in some of their models.
While an innovated and much approiaciated approach, there are some critical areas of Improvement: 1) Data presented in Figure 3. Localization of the plaques could be clearer, especially the subcellular localization. Perhaps could be complemented by alternate imaging methods that show this key property of the model at higher resolution.
2) Data presented in Figure 4. Potentially move oxidative stress validation to supplemental.
3) Data presented in Figure 6. Park et al. tested FDA-approved drugs on both LPLA288T SNP PiB+ iCOs and E4iso iCOs as a validation of their network-based drug screening platform. Results of the validation experiments were not entirely conclusive due in part to lack of a clear dose-dependent relationship. These results therefore would require further validation, potentially by testing with more replicates and in the other cell lines. Drug testing was not performed on PiB+ iCOs that contain the ApoE ε4 allele. In addition to LPLA288T SNP PiB+ iCOs only representing a small subset of the iCOs used in this study as this SNP is not present in other PiB+ iCOs, these iCOs also do not carry the ApoE ε4 allele. Furthermore, although E4iso iCOs were used for drug testing, they are not derived from AD patient iPSCs. Thus, to properly validate their model for precision medicine applications, FDA-approved drugs need to also be tested on PiB+ iCOs that carry the ApoE ε4 allele. These COs better represent the genetic background of human sporadic AD patients than E4iso iCOs, while also carrying the ApoE ε4 allele missing in the LPLA288T SNP PiB+ iCOs.
While I approaciate the unique aspect the mathematical model brings to the paper, I also felt the flow of the paper would be much improved if this was more concisely summarized in the main text.
Reviewer #4: Remarks to the Author: The manuscript submitted describes a high-content drug-screening platform (HCS) that uses human iPSC-derived cerebral organoids (iCOs) in order to identify Alzheimer's disease (AD) drug candidates. For this purpose, pluripotent stem cells (iPSC) derived from participants were used that had been selected on the basis of the presence or absence of pathological features of AD or the ApoE4 genotype. In addition to the use of iCO's from participants with preclinical sporadic AD, a molecular regulatory network model for AD was developed and integrated in order to filter suitable candidates from a library of FDA-approved drugs, the efficacy of which was then tested in the HCS-system. From a technical point of view the described platform contains a couple of interesting and important new features that might help to overcome some of the present limitations in the use of organoids for high-content drug screening purposes. One of them is a stringent emphasis on quality control with respect to producing a large number of homogenous well shaped and evenly sized iCO's. Demonstrating the feasibility of obtaining such iCOs will be relevant for many other applications of organoids and of interest to others in the community.
The selection and comparison of iCO's from participants with or without a high burden of Aβ, total tau and phosphorylated tau deposits, is also an important element of the screening platform presented, providing a potentially more relevant disease model for sporadic AD than those available so far.
Furthermore, the integration of a molecular regulatory network model for AD into the screening platform is put forward by the authors as another key innovative element that will facilitate the identification of suitable drug candidates. The approach nicely exploits multiple sources of information: It takes both gene expression data from iPSC models and prior knowledge from curated AD pathways into consideration, and also prefilters the candidate drugs by removing those that do not pass the blood-brain barrier (BBB).
However, it is not clear that the network analysis in this approach provides much additional and robust filtering information for the candidate drug selection: Instead of performing an unbiased genome-scale network analysis, the authors have already pre-selected a network of only 77 nodes with well-known functional AD associations, mainly from existing curated AD pathways in KEGG AlzPathway. This limits the selection of candidate drug targets significantly -the 77 nodes cover already well-studied AD protein targets, and due to this pre-selection, a subset of them will always be chosen as candidate targets, independent of the observed alterations in the iPSC transcriptomics data.
The known drug associations of these 77 nodes, and the drug BBB-filter will further reduce the number of candidate targets and compounds, again independent of the content of the iPSC data. Thus, before applying any network analysis, the drugs are already pre-filtered to those that have previously been described to target the known AD-associated proteins in the AlzPathway map or the KEGG Alzheimer pathway and that pass the BBB. These will all be reasonable candidate drugs, but with limited novelty.
The network analysis is used as an additional filter to further reduce the candidate drug selection, but it is unclear whether it provides significant and robust filtering information beyond the previous generic filters. The two main reasons why there are doubts regarding the robustness and reliability of this additional filtering step are outlined in the comments 1) and 2) below: The small number of RNA-seq experiments limits the statistical power, the transcriptomics data comes from only 11 participants, which may lead to subject-specific idiosyncasies (1), and the manuscript suggests that the significance scores to determine the differentially expressed genes may not have not been adjusted for multiple hypothesis testing and/or a relaxed p-value cut-off was used.
Main comments: 1.) One limitation is that the data comes from only 11 participants. The data for these subjects may not be representative for the overall population of sporadic AD patients, and subject-specific idiosyncrasies in the data that are not disease-relevant could lead to errors in the transcriptomebased network model and perturbation analysis. The authors mention that they aim at a personalized precision medicine approach, but the distinction between disease-associated biological variance and other sources of variance in the data is still an issue with small numbers of subjects, and small numbers of RNA-seq experiments limiting the statistical power, in particular for the transcriptome-derived network analysis. It is therefore recommended to compare the data against larger public iPSC data for AD, e.g. derived from the GEO database for the transcriptomic analyses.
2.) It is not clear whether the p-values that were used to define the differentially expressed genes (DEGs) were adjusted for multiple hypothesis testing. Using a relaxed p-value cut-off for the pathway analysis (-log10(p-value) > 1, corresponding to a p-value cut-off of 0.1) can be justified by taking into consideration that pathways with an enrichment of many small, close-to-significant alterations are likely disease-relevant. However, since the authors do not only aim at ranking GO processes, but also at identifying individual drug target genes, a false-discovery-rate (FDR) above 10% for individual genes would lead to too many errors (in particular, if the authors did not use FDR-adjusted p-values, but nominal p-values). Therefore, FDR corrected p-values should be used for all analyses that rely on the significance of individual gene alterations.
3) The authors state that the pre-selected network consists of 77 nodes and 203 links. While these nodes representing mainly genes/proteins from KEGG and AlzPathway definitely play important roles in AD, a significantly larger number of genes/proteins will likely be relevant for AD than this pre-filtered subset, and the restriction to mostly KEGG/AlzPathway-derived nodes may bias the results towards target genes in these pathways that are already well-known and whose associated drugs therefore have limited novelty. A possibility to avoid this limitation is to use a pathwayagnostic network analysis (e.g. using a genome-scale network from the STRING web-service or other public resources for genome-scale gene regulatory or protein interaction networks) to identify network clusters of transcriptomic alterations, which are not already captured by the known pathway definitions. 4) To show that the mathematical modeling / network analysis provides a significant added value beyond the compound filtering obtained from the network model pre-selection of 77 nodes and the BBB-filter, it would be useful to compare the ranked target and compound lists with and without the additional network analysis (e.g. testing whether there is an improved enrichment of known AD protein drug targets, such as BACE1, MAOB, MAPT etc., that have been considered in AD clinical trials, in the network analysis derived ranking list).
In summary, the current manuscript describes the use of a drug screening platform that overcomes some of the limitations of current iPSC/iCO human disease models and is of potential interest beyond AD. However, it does not provide sufficient information to show that network analysis improves the filtering significantly beyond what is already achieved by the prior generic filtering steps. Some possibilities to adjust the analysis workflow (e.g. by using a genome-scale network analysis approach) are suggested. in accordance with the BoolNet format, and they have been converted to SBML-qual 22 format using BioLQM v0.6.1. As soon as our paper is accepted, we will add detailed 23 annotations to our model and upload its SBML file to the existing Boolean model 24 repositories such as BioModels, GINsim.org, CellCollective by referencing our study. 25 26

[COMMENT #2] 27
Related to reproducibility, I don't believe the methods for the modeling component are 28 comprehensive enough. For example, what initial state(s) did the authors use for the 29 quality control input-output analyses ( Fig. 5 and Fig. 5-1

), and the actual predictions? 30
Did the authors perform initial condition sensitivity analysis? 31 For the quality control input-output analyses (Supplementary Fig. 7-1, 7-2), the initial 33 state of all but three nodes were set to zero: the value of 'ROS' node was set according 34 to input conditions and the values of 'APOE4' and 'LPL' were fixed depending on the 35 given genetic conditions. For given genetic conditions, qualitative input-output analysis 36 was performed for all 1,000 input levels of 'ROS' between 0% and 100% over 0.1% 37 increment. For each input level, the Boolean network simulation was conducted over 38 1,000-time steps. Each node's state value during the last 300-time steps were tracked, 39 and the average node state value in this period was taken as the 'average node activity'. 40 For the actual prediction (Fig. 5), we used the 'getAttractors' function in the 41 'BoolNet' package. Due to the high computational cost of simulating all initial states, 42 we randomly sampled one million initial states. We confirmed that the main results 43 were not sensitive to the randomly selected initial states through repeated sampling 44 processes. The function's 'method', 'type' parameters were set to 'random' and 45 'synchronous', respectively to find attractors. We also used 'genesON' and 'genesOFF' 46 parameters of the function to simulate genetic conditions and perturbations. 47 Following the reviewer's comment, we have further performed initial condition 48 sensitivity analysis of the input-output analyses. From the total of 77 nodes, a certain 49 number of nodes except 'ROS', 'APOE4' and 'LPL' have been randomly selected and 50 given an initial state value of '1' and the others have been set to '0'. The number of 51 selected nodes has been set from 10 to 70 with 10 nodes interval, and the input-output 52 analysis has been performed for these randomly generated initial state values. The  Table 2). 81 3 82

[COMMENT #4] 83
It is not clear how oxidative stress ("ROS") in the model was used for input-output 84 analyses with varying activity levels, because the model (Fig. 4

and Supplementary 85
Table 2) shows that several other model components regulate ROS --in other words, it 86 is not an independent variable (input) that could be easily used for input-output analyses. 87

[RESPONSE] 88
In our study, oxidative stress ('ROS') was used as an input to describe the effects of 89 aging, which is a major risk factor of Alzheimer's disease that increases pathological 90 protein levels (i.e. beta-amyloid, p-tau, etc.) 2 . However, because biological evidence 91 that ROS is produced by beta-amyloid or Ca 2+ ion exists, the regulatory interactions 92 directed to ROS could not be ruled out in our model 3 . Thus, as the reviewer pointed out, 93 ROS is not strictly an independent variable in our model. Despite this issue, we used 94 ROS as an independent variable to analyze the changes of output node states for fixed 95 initial ROS levels only for input-output analysis. Following the reviewer's comment, we have cited the latest review papers on 117 Boolean/logical modeling in the revised manuscript (line 618-619). A citation for the 118 BoolNet package has also been added (line 627-631). Finally, a citation for the Helikar 119 et al., 2008 has also been added to the 'Methods' section of the revised manuscript, so 120 that appropriate citation for the previous study is made (line 648).

122
Some of the figures are very complex, making them hard to read without zooming in 124 significantly (Fig. 6 is the least readable) 125

[RESPONSE] 126
We have increased the font size on our figures and especially modified Fig. 6. The 127 graph type has been changed in order that the medicinal effects can be seen well at a 128 glance (all of the cell viability data have moved to Supplementary Fig. 11. Instead, 129 PiB + #5 iCOs data were added in the Fig. 6  How did the authors select the drug concentrations? 151

[RESPONSE] 152
Since our drug candidates were subjected to drug repositioning, they had mainly been 153 tested for applications in fields other than neurodegenerative disorders. Most of the 154 references used cell types different from neurons. We could hardly find any studies 155 treating the drugs on brain organoids as it is a relatively novel and unconventional 156 method in drug screening. Consequently, we referred to studies using the drugs in in 157 vitro application and optimized for our system by testing the dosage range that includes 158 the optimal level suggested by the references. 159 When setting the dosage range, we had two points into consideration. Firstly, a lower 160 dosage limit should be tested to observe the potential neurotoxic effects of the drugs. 161 Responsiveness to drugs might differ according to cell types, and iPSC-derived neurons, 162 which resemble primary cells, could be more vulnerable to a higher concentration than 163 conventional cell lines. Secondly, the fact that brain organoids are three-dimensional 164 aggregates of interconnected neurons with complex cytoarchitecture and extracellular 165 matrix should be considered. These might serve as barriers to simple diffusion of drugs 166 and induce essentially different drug penetration dynamics and subsequent uptake by 167 cells compared to conventional monolayer culture. It led us to test the upper 168 concentration limit as equal to or higher than the references. 169 Although the 'Veh' groups (used as controls whereas the treatment groups were used as 187 cases) were not treated with the drugs, they were simultaneously grown in the opti-188 MEM media, which is not suitable for the permanent cell culture of iCOs, during the 189 treatment of drugs with other iCOs. We thought the best time point will be that of 190 'showing a natural tendency of cell death but not too much'. We speculated 48 or 72h 191 (>24h) may be too long for the drug test because the 'Veh' groups also can show 192 excessive loss of neuronal cells (might be called as 'cell senescence') even without the 193 drugs. As expected, the Veh group showed statistically significant loss of their own cell 194 viability from 48h timepoint (reference figure below). So, we decided to use 24h 195 timepoint for the drug treatment. * Line 45: remove "the" before precision medicine/ * Line 50: remove "the" before 200 symptoms/ * Line 76: "in vitro" should be italicized/ * Line 99: "Mathematic" should be 201 "Mathematical"/ * Line 128: Missing "genes" from "differentially expressed genes 202 (DEG) patterns"/ * Line 197: Remove "comprehensive" --while the model is not a "toy 203 model," comprehensive is a relative term, and I don't believe this model is 204 comprehensive (for example many other pathways known to interact with the modeled 205 pathways could be included). Perhaps "relevant" is a more appropriate word/ * Line 313: 206 Similar to above, I recommend removing "the whole" because the model does not 207 include the entire interaction network /* Line 321: "was focusing" does not seem 208 grammatically correct; perhaps change to "was to focus on."?/ * Line 333: remove "but" 209

[RESPONSE] 210
We corrected all of them. Thanks for the kind and detailed comments. However, authors should clearly describe their methods in construction of AD 223 molecular network. [1] How did they narrow genes/proteins/phenotypes down to only 77 224 genes/proteins/phenotypes? It appears to be an arbitrary choice. At least, AlzPathway 225 consists of 1,347 species (genes and proteins) and 129 phenotypes. [2] Authors need to 226 explain the reason why they focus on MAPK signaling pathway, WNT signaling 227 pathway, and PI3K-AKT signaling pathway as molecular regulatory network of AD 228 (Fig. 4A) which is a basis for network analysis in this manuscript. [3] By the way, Fig.  229 4A is the correct figure? For example, Reelin has relationships with not only Apoer2 but 230 also VLDLR and Apoptosis according to a logic table (Supplementary Table 2). 231 [4] Authors also should clearly describe their methods in analysis of the AD network 232 model and identification of candidate drugs. For example, Fig. 5A illustrates up-233 regulated and down-regulated pathways according to their perturbation analysis, but 234 they did not explain the definition of "up-regulated" and "down-regulated" 235 pathways, and [5] they did not show their results of perturbation analysis. [6] Authors also 236 should clearly describe their methods of attractor landscape analysis and unfortunately 237 [7] attractor landscape drawn in Figure 5B is too small to see. [8] By the way, the idea of 238 phenotype score is interesting to estimate the pathological level, but there are several 239 concerns. Authors choose key proteins and phenotypes such as Aβ, p-tau, synapse loss, 240 apoptosis, and autophagy for calculation of phenotype score, but how to choose these 241 key proteins and phenotypes? [9] Phenotype score looks work well but p-tau proteins 242 have double impact on the score because p-tau activates synapse loss. They developed a 243 network model but they did not consider these kinds of network effects in the 244 calculation of phenotype score. If phenotype score well works to estimate the 245 pathological level, [10] why did they need to develop a network model and conduct a 246 perturbation analysis? Why don't authors directly conduct the phenotype score 247

calculation? 248
Basic idea of this manuscript is a novel and interesting in developing drug-screening for 249 Alzheimer's disease, but they need to clearly describe their methods and results. 250 251

[COMMENT #1] 252
How did they narrow genes/proteins/phenotypes down to only 77 253 genes/proteins/phenotypes? It appears to be an arbitrary choice. At least, AlzPathway 254 consists of 1,347 species (genes and proteins) and 129 phenotypes. 255

[RESPONSE] 256
The AlzPathway has various types of components (e.g. gene, protein, signaling pathway, 257 etc.), so its components could not be directly used as nodes for the network model. For 258 this reason, we used the AlzPathway for sorting out the list of signaling pathways. Of 259 the components in the AlzPathway, only the signaling pathways of the neuron 260 (indicated by 'n_' prefix) were considered for the construction of the network model. 261 These pathways were grouped into Notch, RELN, MAPK, Jak/Stat, Wnt, NR1/NR2R, 262 Ca 2+ , TGFβ and mTOR (PI3K-AKT) signaling pathways. Among them, TGFβ signaling 263 pathway related to microglia 8 was excluded. Based on these signaling pathway 264 information, we searched for experimental studies in neuronal context to the best of our 265 knowledge. By using this literature-based approach, we reconstructed simple regulatory 266 links (i.e. activate or inactivate) and actual regulatory relationships such as update logic 267 of node state (e.g. A node = B node AND C node) (Supplementary Table 2). 268 269

[COMMENT #2] 270
Authors need to explain the reason why they focus on MAPK signaling pathway, WNT 271 signaling pathway, and PI3K-AKT signaling pathway as molecular regulatory network 272 of AD (Fig. 4a) which is a basis for network analysis in this manuscript. 273

[RESPONSE] 274
In our study, we focused on APOE4, LPL-related signaling pathways in line with the 275 experimental results in the main manuscript (Fig. 2). Through a literature-based 276 approach, we constructed the model network primarily based on MAPK, Wnt, and 277 mTOR signaling pathways 9-11 , and also considered other signaling pathways in the 278 AlzPathway. We have added more explanation on the model construction in the revised 279 manuscript (line 608-612). 280 281

[COMMENT #3] 282
By the way, Fig. 4a is the correct figure? For example, Reelin has relationships with not 283 only Apoer2 but also VLDLR and Apoptosis according to a logic table 284 Table 2). 285

[RESPONSE] 286
As the reviewer pointed out, there was a typographical error in Authors also should clearly describe their methods in analysis of the AD network model 292 and identification of candidate drugs. For example, Fig. 5a illustrates up-regulated and 293 9 down-regulated pathways according to their perturbation analysis, but they did not 294 explain the definition of "up-regulated" and "down-regulated" pathways. 295 Authors also should clearly describe their methods of attractor landscape analysis 296

[RESPONSE] 297
As we described in the 'Mathematical modeling' section of 'Methods', attractor is a 298 network state that represents a specific biological phenotype. Attractor landscape is the 299 landscape of all attractors and the set of all states converging to each attractor (basin of 300 attraction) in the state space of a given network model. Attractor landscape analysis 301 enables quantitative evaluation of the network system by transforming complex 302 dynamical properties of the network model into the states of convergence (attractors) 303 and the converging propensity of initial states (basins). Using attractor landscape 304 information, we can quantify node activities by averaging attractor states weighted by 305 their basin sizes. In other words, attractor landscape analysis is a kind of procedure that 306 can quantify node activities using attractor landscape information and further estimate 307 the perturbation effect in silico. We have added the aforementioned description to 308 'Attractor landscape analysis' in 'Methods' of the revised manuscript (line 633-641).

310
The term 'up-regulated' and 'down-regulated' means the increase and decrease, 311 respectively, of representative nodes (genes or proteins) for each signaling pathway. For 312 instance, 'up-regulated' and 'down-regulated' correspond to the increase and decrease, 313 respectively, of ERK node activity in case of the MAPK signaling pathway. The 314 representative nodes for each signaling pathway are described on the right-side table of 315 the Supplementary Fig. 8. 316 317

[COMMENT #5] 318
They did not show their results of perturbation analysis. 319

[RESPONSE] 320
The major perturbation analysis results used in the experiments (Fig. 5)   In the conceptual image of the '① Node perturbation' box of Fig. 5a, each grey circle 332 in 'Attractor landscape' represents network state, which is a collection of node states (1 333 or 0). The black circle is the attractor to which all states within the basin of attraction 334 converge. The grey circles within the boundary of attractor A belong to the basin of A. 335 The basin size of attractor A, which is the number of grey circles in the basin of A, is 336 larger than the basin size of attractor B. 337 338

[COMMENT #8] 339
By the way, the idea of phenotype score is interesting to estimate the pathological level, 340 but there are several concerns. Authors choose key proteins and phenotypes such as Aβ, 341 p-tau, synapse loss, apoptosis, and autophagy for calculation of phenotype score, but 342 how to choose these key proteins and phenotypes? 343 The reasons why we chose Aβ, p-tau, neuron loss, synapse loss, and autophagy as key 345 proteins and phenotypes are as follows. First, because Aβ and p-tau are the well-known 346 biological markers of Alzheimer's disease, they were chosen as the target proteins the 347 expression levels of which should be decreased. Second, the goal of our study is to 348 identify an optimal target that can decrease the neuron/synapse loss, which is associated 349 with the cognitive impairment in Alzheimer's disease, by analyzing the network model. 350 Third, based on previous findings that the degradation of autophagy function increases 351 accumulation of Aβ, we included autophagy as a key phenotype 12,13 . Choosing 352 autophagy as a key phenotype helps to understand whether the accumulation 353 mechanism of Aβ is due to an increase of Aβ production or a decrease in autophagy 354 function. For these reasons, we chose Aβ, p-tau, neuron loss, synapse loss, and 355 autophagy as key proteins and phenotypes in this study. 356 As we described in the 'Analysis of the AD network model and identification of 357 candidate drugs' part of 'Results', we calculated the phenotype score using the key 358 proteins and phenotypes' node activities obtained from attractor landscape analysis of 359 our network model. Therefore, simulation of the network model is a prerequisite to 360 obtain the phenotype score. In the network model with a given genetic condition (e.g. 361 APOE4, LPL), we can identify which node contributes to producing a desirable 362 phenotype score using the node perturbation analysis as described in 'Result' of the 363 main manuscript. The desirable phenotype score refers to the case where the activities 364 of the pathological proteins and phenotypes are low. Hence, both perturbation analysis 365 and calculation of the phenotype score are required to find the optimal target(s) for a 366 given genetic condition. 367 368

[COMMENT #9] 369
Phenotype score looks work well but p-tau proteins have double impact on the score 370 because p-tau activates synapse loss. They developed a network model but they did not 371 consider these kinds of network effects in the calculation of phenotype score. 372

[RESPONSE] 373
Although there is a positive link from 'p-tau' to 'Synapse loss', activation of 'p-tau' 374 alone cannot satisfy the sufficient condition for activating 'Synapse loss' according to 375 the authors do not only aim at ranking GO processes, but also at identifying individual 501 drug target genes, a false-discovery-rate (FDR) above 10% for individual genes would 502 lead to too many errors (in particular, if the authors did not use FDR-adjusted p-values, 503 but nominal p-values). Therefore, FDR corrected p-values should be used for all 504 analyses that rely on the significance of individual gene alterations. 505

[RESPONSE] 506
Thank you for the critical comment. It needs to be clarified that different correction 507 strategies were adopted for each RNA-seq analysis. 508 First, we already applied the FDR corrected p-values for the validation of mathematical 509 model and identification of individual drug target genes ( Fig. 5b and Supplementary  510   Fig. 8,9). We have mentioned this in the legend of Supplementary Fig. 8. 511 Second, to demonstrate that gene expression profiles (GO analysis) differ between PiB -512 and PiB + groups (Fig 2i and Supplementary Fig 2a, c, f), we subjected five samples 513 from each group, performed RNA-seq on each sample and compared gene expression 514 profiles (GO analysis) by using the selected DEGs. For this specific analysis, we did not 515 apply multiple correction for the selection of DEGs because we intended to be more 516 inclusive in documenting subtle differences in gene expression changes, as the reviewer 517 pointed out. However, since we had not applied FDR correction, even for the GO 518 analysis deriving significant GO terms from the selected DEGs, we have changed p-519 values to FDR-adjusted p-values (Fig 2i and Supplementary Fig. 3). We have 520 mentioned this in the legend of Fig 2 and Supplementary Fig. 3  521 Regarding the reviewer's comment, we have also clearly mentioned it in Methods 522 section (line 592-595, line 600-603), please refer to the Methods section. 523 524

[COMMENT #3] 525
The authors state that the pre-selected network consists of 77 nodes and 203 links. 526 While these nodes representing mainly genes/proteins from KEGG and AlzPathway 527 definitely play important roles in AD, a significantly larger number of genes/proteins 528 will likely be relevant for AD than this pre-filtered subset, and the restriction to mostly 529 KEGG/AlzPathway-derived nodes may bias the results towards target genes in these 530 pathways that are already well-known and whose associated drugs therefore have 531 limited novelty. A possibility to avoid this limitation is to use a pathway-agnostic 532 network analysis (e.g. using a genome-scale network from the STRING web-service or 533 other public resources for genome-scale gene regulatory or protein interaction networks) 534 to identify network clusters of transcriptomic alterations, which are not already captured 535 by the known pathway definitions. 536

[RESPONSE] 537
The AlzPathway has various types of components (e.g. gene, protein, signaling pathway, 538 etc.), so its components could not be directly used as nodes for the network model. For 539 this reason, we used the AlzPathway for sorting out the relevant signaling pathway lists. 540 Of the components in the AlzPathway, only the signaling pathways of the neuron 541 (indicated by 'n_' prefix) were considered for the construction of the network model. 542 These pathways were grouped into Notch, RELN, MAPK, Jak/Stat, Wnt, NR1/NR2R, 543 Ca 2+ , TGFβ and mTOR (PI3K-AKT) signaling pathways. Among them, TGFβ signaling 544 pathway related to microglia was excluded. 545 Based on these signaling pathway information, we searched for experimental studies 546 only in neuronal context to the best of our knowledge. By using this literature-based 547 approach, we reconstructed simple regulatory links (i.e. activate or inactivate) and 548 actual regulatory relationships such as update logic of node state (e.g. A node = B node 549 AND C node) (Supplementary Table 2). 550 In our study, we focused on APOE4, LPL-related signaling pathways in line with the 551 experimental results in the main manuscript (Fig. 2). Through a literature-based 552 approach, we constructed the model network primarily based on MAPK, Wnt, and 553 mTOR signaling pathways 9-11 , and also considered other signaling pathways in the 554 AlzPathway. We have added more explanation to the model construction in the revised 555 manuscript (line 608-612). 556 The purposes of our study are as follows: 1) to understand complex dynamical 557 properties of the regulatory network, which cannot be done by simple correlation 558 analysis among the network components, and 2) to identify optimal candidate targets for 559 each genetic condition based on quantitative dynamical analysis of the regulatory 560 network model. As the reviewer pointed out, pathway-agnostic network analysis may 561 have benefits in finding novel targets. However, since we aim to find the optimal targets 562 for a given genetic condition rather than to find novel targets, the proposed approach 563 does not fit our research objectives. Also, public resources such as STRING contain 564 collective information from various cell-types, so we considered that such kinds of 565 resources might not be suitable for constructing our neuron-specific network model. 566 Nevertheless, the reviewer's comment will be very helpful in our future study to find 567 out novel targets by further extending our network model to include diverse genetic risk 568 factors and components. We appreciate the reviewer's valuable comment in this respect. 569 570

[COMMENT #4] 571
To show that the mathematical modeling / network analysis provides a significant added 572 value beyond the compound filtering obtained from the network model pre-selection of 573 77 nodes and the BBB-filter, it would be useful to compare the ranked target and 574 compound lists with and without the additional network analysis (e.g. testing whether 575 there is an improved enrichment of known AD protein drug targets, such as BACE1, 576 MAOB, MAPT etc., that have been considered in AD clinical trials, in the network 577 analysis derived ranking list). 578

[RESPONSE] 579
There were several network analysis studies that constructed gene or protein interaction 580 networks using databases such as the STRING and found network clusters/modules of 581 transcriptomic alterations that are significantly related to diseases. However, the data 582