Physics-informed deep learning characterizes morphodynamics of Asian soybean rust disease

Medicines and agricultural biocides are often discovered using large phenotypic screens across hundreds of compounds, where visible effects of whole organisms are compared to gauge efficacy and possible modes of action. However, such analysis is often limited to human-defined and static features. Here, we introduce a novel framework that can characterize shape changes (morphodynamics) for cell-drug interactions directly from images, and use it to interpret perturbed development of Phakopsora pachyrhizi, the Asian soybean rust crop pathogen. We describe population development over a 2D space of shapes (morphospace) using two models with condition-dependent parameters: a top-down Fokker-Planck model of diffusive development over Waddington-type landscapes, and a bottom-up model of tip growth. We discover a variety of landscapes, describing phenotype transitions during growth, and identify possible perturbations in the tip growth machinery that cause this variation. This demonstrates a widely-applicable integration of unsupervised learning and biophysical modeling.

Supplementary Figure 1: Convergence of the physics-informed neural network (PINN). (a) An ablation analysis comparing how L total (running mean over 200 mini-batches) decreases as training progresses reveals that a diffusivity with both spatial and time dependence is the best model. The majority of the benefit likely comes from the dynamics in the spore region of morphospace, where diffusion is very high at first, and then strongly decreases such that not all spores germinate. (b) We repeat PINN network training three times for each condition, with different mini-batches, and L total (running mean over 200 mini-batches) is shown here for each repeat of DMSO. (c) An approximation of the dynamics of the fraction of spores (≈ sporep (x, t)) across snapshots for DMSO (with the first snapshot not shown due to its much higher fraction of spores), found by numerically integrating a box around the spore PDF peak. The data is shown in red, and the three repeats have the same coloring as in (b). The PINN first explores smooth low-frequency solutions, fitting trends common to all snapshots, before ultimately beginning to overfit to the individual snapshots, as shown for one repeat at 60 hours in black. We stop training when the PINN begins to fit to the individual snapshots, which approximately corresponds to 30, 30, 30, 20, 25 and 25 hours for DMSO and Compounds A, B, C (0.041 mgL -1 ), C (10 mgL -1 ) and X, respectively. (d) Sketch of the loss landscape, whereby the global minimum is an overfitted solution, and there may be many equally good solutions before overfitting. (e) The landscapes from each of the three repeats after 30 hours of training for DMSO show many common features in the central data-rich region. (f ) For each of the conditions with significant germination (i.e. excluding Compound C at 10 mgL -1 ), three outputs are shown: the uncertainty in the force magnitude, F , calculated from the standard deviation across the three training repeats; the mean σ (from Eq. 1), averaged over time for the same training repeat as those of the landscapes shown in the other figures, and only calculated over regions where the PDFs are above 10 −3 ; and the uncertainty in σ, calculated in the same way as the uncertainty for F . All outputs are expressed in terms of morphospace units, m.u. Source data are provided for (c). The same as described above, but for Compound C at 0.041 mgL -1 , alongside the landscapes for Compound C at 0.041 mgL -1 (with contours along equal landscape values, spaced 0.14 m.u. 2 min -1 apart, where m.u. stands for morphospace units) and 10 mgL -1 (with contours spaced 0.02 m.u. 2 min -1 apart). Supplementary Figure 3: Validation of the landscapes and diffusivities learned by the PINN. For each condition (a-f ), three panels are shown: the first panel is the data kernel density estimate (KDE), the second is the KDE over simulations, and the third is the error (data KDE -simulation KDE). All are displayed on a logarithmic scale, and the error is truncated at 10 −2.5 , which is the probability density generated by a single particle, in order to highlight more systematic errors. For the forward simulations, particle starting positions were sampled from the initial probability distribution learned by the PINN, and then simulations were run by evaluating the potential and diffusivity on a 1000 × 1000 spatial grid, with 20 snapshots in time for the diffusivity. The figure shows good agreement across all conditions, validating the landscapes and diffusivities learned by the PINN. Supplementary Figure 4: Comparison of simulation and data trajectories. (a) An example sequence from the DMSO time-lapse videos. Images were taken every 3 min, from 60 min after mixing with the compounds. (b) Trajectories of sequential frames of the time-lapse videos (colored) and a sample of DMSO simulations (black). For the forward simulations, particle starting positions were sampled from the initial probability distribution learned by the PINN, and then simulations were run by evaluating the potential and diffusivity on a 1000 × 1000 spatial grid, with 20 snapshots in time for the diffusivity. (c) Mean squared displacement (MSD, in terms of morphospace units, m.u.) plots against time for the time-lapse videos (colored) and forward simulations (black). Time-lapse videos were taken under higher temperatures, which results in early germination. A confusion matrix of the mean absolute differences (d) of the plots is also shown. For each simulation-data pairing (videos down the rows, simulations across columns), the time series were shifted horizontally and the result for each pairing taken to be the minimum of the mean absolute differences across the shifts. Simulations match their corresponding data generally the best, except for the simulations of Compound X. (d) The entropy of PDFs from a KDE over single particle simulations of Eq. 1 reveals that entropy always increases with time. Entropy is calculated as − x1 x2 p(x) log 2 p(x)∆x 1 ∆x 2 with x = (x 1 , x 2 ) and only summing over morphospace regions where the PDFs are above 10 −3 . Source data are provided for (c-d).

Supplementary Figure 3: Correspondence between the landscapes and morphospace
Supplementary Figure 5: Data-driven development of the tip growth model  Variation in global direction, θ global , for the three tip bending models tested (using bending MAP values for Compound A, with a growth rate of 0.75 µm min -1 ). Model 1 is a random walk in θ global , Model 2 is a random walk in path curvature, κ, and Model 3 is a persistent random walk in κ, with relaxation to straight growth. (c) Model selection using ABC-SMC. For early populations where the acceptance threshold is high, the lower dimensional parameter spaces of models 1 & 2 lead to better fits. At lower acceptance thresholds, however, models 2 & 3 fit better, validating model conception in the tip frame, and ultimately the relaxation to straight growth in Model 3 is required to reproduce the data distribution. (d) Comparison of Compound A snapshot data (pink) and MAP simulations (black) for the three models at 210 min, with an enlarged example of a randomly selected simulation and data fungus shown in the inset of the Model 1 box. While all models introduce bending too early for some fungi (the region below the spore where there are simulations but no data), Model 3 can reproduce the feature distribution best. In particular, it is the only model that can reproduce fungi with multiple bends in alternating directions. Overall, this feature is less well separated in this 2D morphospace that prioritizes global features. Source data are provided for (a, c).

31
After determining the spore concentration by using a hemocytometer (Neubauer improved), the required 32 volume of two-fold concentrated suspension with 20,000 spores per mL in 0.0015 % Tween 20 was prepared.

33
The assay was started immediately by mixing 50 µL of spore suspension with treatment solution in the 96-well 34 plates prepared in advance, resulting in a final one-fold concentration of the test compounds (Carbendazim:  (NanoEnTek Inc.). The germination assay was set up as described above, but after mixing spores with 52 treatment solutions the assay plates were incubated directly in the imaging device, which itself was placed 53 inside of a climate cabinet set to 20°C. Imaging started with an autofocus run, followed by 321 runs with 54 3 min intervals to cover germination between +1 and +17 h after setup. No replicate wells were prepared 55 to minimize interval times, but 5 different spots were imaged per well. The transmission light LED and the 56 10x objective were used with settings: 15 ms exposure time, LED power "3", brightness correction "10".

57
Image Processing

58
The size of the snapshot image sets necessitated fully-automated processing. We extracted fungus contours 59 using adaptive binarization, where the threshold value varies based on the statistics of a surrounding window. 60 We used 100×100 windows (the approximate size of lighting defects) out of full images of 503×685 pixels, and in opposite directions until the spore was hit. We replaced all spores with identical circles, so as to prioritize 82 modeling of the germ tube; the resulting morphospace point is then widened into a spore region through the 83 kernel density estimation. 84 We removed overlapping fungi, which we note means larger morphologies are more likely to be removed  For the time-lapse videos, we again binarized (using trial and error for a suitable threshold value, and 94 without adaptive binarization, since there were not significant lighting defects), and found the fungus con-95 tours. We then found series of contours across frames whose center of mass were closest, and manually looked 96 through these to find those that corresponded to tracking an individual fungus. We then used ImageJ to  For the PINN, the loss function to be minimized comprizes four terms, with the first three calculated over 132 random mini-batches of N data points, and the final one over the full spatial grid of M data points. The 133 first is the mean squared difference between the learned PDF,p(x j , t j ), and data, p(x j , t j ), with {x j , t j } in the nine snapshots. The second is the mean squared PDF at the boundary, with {x j , t j } selected from 10 6 uniformly distributed boundary points, and the third term is the mean 136 squared PDE residual (N , given in Eq. 3), with {x j , t j } selected from 10 6 points uniformly distributed over the whole domain. The final term ensures 138 the PDF integrates to one: with x j covering the full spatial grid and t randomly selected.
For the total loss (Eq. 4), we used hyperparameters of 1, 1, 500, 0.01 for a, b, c and d. Previous work 141 using PINNs to solve the Fokker-Planck equation [6,7] used a weighting for L PDE of 100, and the same 142 values for the other hyperparameters as we used. Since our data was of snapshots of different spore batches, 143 we increased the weighting to 500, to prioritize more the PDE fitting over the data. While initially selected 144 by trial and error through inspecting single particle forward simulations, this choice was later vindicated only on the single-snapshot level, but also across many snapshots; for example if two subsequent snapshots 165 have similar low-probability latent variables, this will be captured by the PINN. However, such patterns 166 become increasingly less problematic as the number of snapshots increases. An alternative solution is to 167 constrain the solution physically, as we did with the tip growth model.

168
The three neural networks had 5 fully connected layers, each with 50 neurons, with residual skip connec-169 tions, and swish activations between layers. A softplus output activation was used for the PDF, and sigmoid 170 was used for the potential and diffusivity, with the potential multiplied by 3 to give a [0, 3] output range.

171
This sigmoid constraint prevents unphysical solutions that can arise with unbounded force and diffusivity.

172
Output variables that share inputs (e.g. the PDF and diffusivity) can be outputted from a single neural 173 network if they are likely to have similar features, for increased computational efficiency. We used the Adam The minimal model for fungus growth was composed of two equations: one for lengthening and another for 198 tip bending, and we ran model selection on three candidate models for the bending part.

199
For fungus lengthening, we used the 3-parameter lognormal distribution to model both germination time, t g , and growth rate, α. The probability density function of the 3-parameter lognormal distribution is given where σ is a shape parameter, s is a scale parameter (also the median), and loc is a location parameter (the 203 lower bound). The 2-parameter distribution has loc set to zero. 204 We modeled germination time, t g , as distributed according to t g ∼ lognormal(s tg , σ tg , loc tg ), and growth 205 rate, α, as distributed according to α = loc α − x, with x ∼ lognormal(s α , σ α , 0) and resampling for negative  We compared three models for tip bending, and used ABC-SMC for both model selection and parameter 217 inference. As described in the main text, the comparison was done using both the fungus length data and 210 218 min morphospace embeddings, and Model 3 was found to reproduce the data best. For all of the following, 219 σ is a noise parameter that was fitted, and dW is the Wiener process. Model 1 was a random walk in the 220 global direction, θ global , a simple model commonly used in the literature: Model 2 was a random walk in the curvature of the growth path, κ, in order to connect to cell tip mechanics: Model 3 was a persistent random walk in the curvature, with an additional parameter for relaxation to 224 straight growth, τ −1 , motivated by work analysing fission yeast tip growth mechanics [9]: Models 2 and 3 can be loosely connected to a diffusing growth zone at the tip as has been described in Coordinates were then converted to images using the polylines function in OpenCV.
For selecting the optimal bending model, we ran ABC-SMC with a population size of 40, each with 1000 235 simulations, for 9 steps. All prior distributions were uniform distributions over the following ranges The compounds were identified at pre-screening by eye to show a range of phenotypes.

241
Compound A (methyl benzimidazol-2-ylcarbamate) is a widely-used fungicide that inhibits the assembly 242 of tubulin subunits into functional microtubules, which are an essential part of the cytoskeleton [10]. Micro-243 tubules participate in maintaining the shape of cells, the distribution of organelles, the transport of materials it inhibits fungal PI3K, then the observed phenotype could be due to the disturbance of cellular signaling 251 required for normal tip growth.

252
Compound C (benzovindiflupyr) has complex II of the respiratory chain as the target, and consequently 253 an inhibition can lead to a depletion of energy. Therefore, depending on the concentration, the cell will stop 254 growth because it lacks the ability to produce the required metabolites.