Computational image analysis reveals the structural complexity of Toxoplasma gondii tissue cysts

Toxoplasma gondii is an obligate intracellular parasite infecting up to one third of the human population. The central event in the pathogenesis of toxoplasmosis is the conversion of tachyzoites into encysted bradyzoites. A novel approach to analyze the structure of in vivo-derived tissue cysts may be the increasingly used computational image analysis. The objective of this study was to quantify the geometrical complexity of T. gondii cysts by morphological, particle, and fractal analysis, as well as to determine if it is impacted by parasite strain, cyst age, and host type. A total of 31 images of T. gondii brain cysts of four type-2 strains (Me49, and local isolates BGD1, BGD14, and BGD26) was analyzed using ImageJ software. The parameters of interest included diameter, circularity, packing density (PD), fractal dimension (FD), and lacunarity. Although cyst diameter varied widely, its negative correlation with PD was observed. Circularity was remarkably close to 1, indicating a perfectly round shape of the cysts. PD and FD did not vary among cysts of different strains, age, and derived from mice of different genetic background. Conversely, lacunarity, which is a measure of heterogeneity, was significantly lower for BGD1 strain vs. all other strains, and higher for Me49 vs. BGD14 and BGD26, but did not differ among Me49 cysts of different age, or those derived from genetically different mice. The results indicate a highly uniform structure and occupancy of the different T. gondii tissue cysts. This study furthers the use of image analysis in describing the structural complexity of T. gondii cyst morphology, and presents the first application of fractal analysis for this purpose. The presented results show that use of a freely available software is a cost-effective approach to advance automated image scoring for T. gondii cysts.


Introduction
Toxoplasma gondii (T. gondii) is an obligate intracellular parasite belonging to the phylum Apicomplexa that can infect all warm-blooded animals, including humans. It is considered as one of the most successful parasites on the Earth, estimated to infect up to one third of the human population [1]. While the infection is generally mild and self-limiting in immunocompetent individuals, it has the potential to cause severe outcomes in immunocompromised patients and the fetus [2].
The three infectious life stages of T. gondii are rapidly proliferating tachyzoites, slowly growing bradyzoites (inside tissue cysts), and sporozoites (inside oocysts). During its complex life cycle, tachyzoites differentiate into bradyzoites, which in an immunocompetent host form tissue cysts, predominantly in the central nervous system and muscles [3]. Encystment is recognized as a central event in the pathogenesis, persistence, and transmission of toxoplasmosis.
Studies focusing on tissue cyst biology have been conducted for decades. Despite numerous in vitro studies [4][5][6], in vivo research that may provide more reliable observations on the structure and organization of T. gondii cysts is still limited [7]. In particular, precise quantitative evaluation of parameters describing cyst structure can help gain additional insights about the cyst structure and development.
In recent years, different types of computational image analysis have found their application in life sciences, with the potential to change the future of medicine [8,9]. Fractal analysis, as a method for quantifying the complexity of natural objects, is being increasingly used in biomedical research [10][11][12][13], including a few studies in the field of microbiology, focusing on bacteria [14,15], fungi [16], and Plasmodium [17]. In addition, particle analysis, a computational biology method for precise cell counting [18], has been successfully applied in several studies for the counting of parasites, including T. gondii [7,19]. However, fractal analysis has not yet been used for studying T. gondii. The ability of fractal analysis to describe natural processes in various tissues implies that it may be used for studying the structural complexity of T.gondii and eventually facilitate research on in vivo-derived tissue cysts.
We here used the freely available ImageJ software (NIH, Bethesda, MD, USA) to quantify the geometrical complexity of T. gondii cysts by fractal analysis and to determine if it is impacted by a number of parasite or host type. Additionally, morphological and particle analyses of T. gondii cysts were applied to gain further insight into the cyst shape uniformity, as well as to a possible correlation between the cyst size and the number of parasites.

Study design
Images of T. gondii cysts analyzed in this paper were obtained in two ways, by selection from the photo archive of strains isolated and maintained in the Serbian National Reference Laboratory for Toxoplasmosis (NRLT), and from cysts harvested from mice infected specifically for this study.
Four different T. gondii strains belonging to genotype-2 were included.

Archived cyst image source
Three out of the four analyzed strains, referred to as BGD1, BGD14, and BGD26, were isolated from human biological materials and thereafter maintained by regular passages in Swiss Webster (SW) mice by intraoesophageal gavage of brain suspensions containing cysts of the respective strains [20][21][22]. In this study, archived images of cysts obtained from the 8 th (BGD1), 21 st (BGD14), and 2 nd (BGD26) passage, were used. These cysts were more mature than those from experiment, obtained from mice infected 6 to 15 months previously.

Cysts obtained from experimental infection
For experimental infections, two genetically different hosts, SW and BALB/c mice (Medical Military Academy Animal Research Facility, Belgrade, Serbia) were used. All mice were female and weighed 18-20 g at the beginning of the experiment. Mice were housed five per cage at the Institute for Medical Research Animal Research Facility, and offered regular mouse feed and drinking water ad libitum. Mice were inoculated intraperitoneally (i.p.) with 10 4 tachyzoites of the reference Me49 strain and the BGD1 strain (from the 32 nd passage). Tachyzoites were harvested from cultivated Vero (ATCC No. CCL-81) cells inoculated with peritoneal exudates of mice, previously infected i.p. with cysts. Me49 parasites were used for infecting both SW and BALB/c mice and cysts were harvested at two time points, at 6 and 12 weeks post-infection (p.i.). BGD1 parasites were inoculated only into SW mice, and cysts were harvested 8 weeks p.i.
To harvest brain cysts, mice were sacrificed by cervical dislocation, brains removed and homogenized with a syringe in 1mL of saline each.

Visualization of T. gondii cysts
For cyst visualization, 25 μl of the brain suspension was placed on slides without fixation and examined under a phase-contrast microscope (Axioskop 2 Plus Zeiss) magnification of 1000x (oil immersion). Cysts were morphologically recognized and photographed by Zeiss Axio-CamMRc 5 in a live image window. Images of both archived and experimentally derived cysts were taken in the same way. Cyst images were saved as 2584 x 1936 px 2 TIFF files. Image format was 24bit RGB. Several images of one cyst were taken, and the best-focused one was selected for further evaluation.
Diameter was measured immediately after taking the images of microscopically visualized cysts.

Image analysis
A standardized analysis protocol using ImageJ was applied on all selected high resolution images, details of which are given in the Supplementary data-S1 Appendix.
Circularity and particle analysis. Images of T. gondii cysts, saved as the grayscale TIFF, were studied by means of morphological, particle, and fractal analysis. All analyses were performed using the ImageJ software package (NIH, Bethesda, MD, USA). The main morphological parameter [23,24] of interest was cyst circularity, defined as Circularity is a measure of how close a shape is to a perfect circle. Shapes close to a circle have circularity values close to 1.0, whereas elongated shapes have very low circularity (close to 0.0). Cysts were carefully manually delineated, as shown in Fig 1A. Parameters including Area, Perimeter, and Circularity were accurately determined by the ImageJ software. The delineated cyst cross-sectional areas were transferred to black background images for subsequent particle analysis.
For particle analysis, images were first enhanced using histogram equalization. During this operation, the distribution of the most frequent pixel intensity values was spread out, resulting in the global contrast enhancement. The resulting image, corresponding to the original one in Fig 1A, is shown in Fig 1B. These enhanced contrast images were subjected to a local autothresholding, using the Bernsen method with a window radius equal to 7 pixels. The Bernsen method belongs to a group of locally adaptive thresholding methods based on contrast. The ImageJ implementation uses circular windows, which is convenient for the analysis of T. gondii cysts. The obtained binary (black and white, B/W) images were further segmented by a watershed algorithm, as shown in Fig 1C, to divide the larger foreground areas of black pixels into the particles to be counted. Segmentation by a watershed algorithm is typically used to enable precise counting of closely packed cells [18,19].
It should be noted that the light areas corresponding to bradyzoites, visible in the cyst cross-section in Fig 1B, were shown as foreground black pixels in Fig 1C. These foreground pixel areas were further segmented by a watershed algorithm for two reasons: a) it seems that, due to a low local border contrast, multiple bradyzoite areas were recognized as a single entity by Bernsen thresholding, and b) the particle counting procedure, currently available in ImageJ, is known to work properly only with clearly separated particle areas. In that manner, all of the cyst cross-section images were subjected to an automated processing and the obtained particle numbers, N p , were expected to be proportional to the real numbers of bradyzoites, with very close proportionality coefficients. The parameter of interest, expected to be proportional to the actual count, was termed the 'packing density' (PD), and calculated as the number of counted particles divided by the cross-sectional area. The corresponding unit is μm -2 .
It is very difficult to claim with certainty where the actual bradyzoite contours are located in an image but the application of the above methods allows for a highly objective comparison of the analyzed specimens. Namely, carrying out the identical procedures in all analyzed images is expected to result in very similar proportionality coefficients with respect to the actual number of particles.
Fractal analysis. The visual texture of overlapping bradyzoites was quantified using the fractal dimension (FD) as a first-order textural complexity parameter and fractal lacunarity as an additional second-order parameter. To ensure that all cysts were analyzed in a consistent manner, cutouts sized 450 pixels by 450 pixels were prepared from the central parts of the cross-sections of all cysts, as shown in Fig 2A. The corresponding areas of the specimens equaled 580.325 μm 2 . Binary images for fractal analysis were obtained by global thresholding using two close threshold values in the vicinity of the Otsu threshold, which is equal to the average of mean grayscale levels of the two sets of pixels being partitioned. FD and lacunarity of the extracted contours of the bradyzoite areas, shown in Fig 2B, were calculated by the boxcounting method and cumulative mass method. Both methods calculated the FD with high accuracy and gave the same results.
The increase in the complexity of a fractal at the higher resolution scale has been shown to follow the power scaling laws. The relationship between a considered measure η(l) and a measurement scale, l, typical for multiplicative processes is known to be of a following form [8,25,26]: Both calculation methods, applied here, rely on producing the logarithm-logarithm plot of a measure versus the scale or grid cell size and fitting the regression line based on the obtained points. The calculation of the exponent, FD, based on the regression line slope is reliable and not too sensitive to the fluctuations of A(l). We generated a series of grid cell sizes using the scaled series 7/8 (seven over eight), meaning that the ratio of consecutive grid cell sizes is about 0.875. The following grid cell sizes were used : {17, 19, 22, 25, 29, 33, 38, 43, 49, 56, 64, 74,  84, 96, 110, 126, 144, 164, 188, 214, 245, 280, 320} pixels. To reduce the calculation variability due to a particular grid cell positioning, twelve randomly positioned grids per cell size per image were used. An example positioning of a grid of cell size 38 pixels is shown in Fig 2B. In the box-counting method, FD is estimated from a change in the number of non-empty grid cells, N(ε)~ε FD , with a change of the cell size normalized with respect to the image perimeter L, ε = l/L, as illustrated in Fig 2C. The FD was obtained as a negative slope of the best-fit regression line: Each FD calculation is accompanied by the correlation coefficient, r 2 , describing a goodness of fit of the regression line, in a particular case. In the cumulative mass method, a probability of finding a certain pixel number within a given cell was estimated by counting the number of pixels contained in each grid cell [25]. In this case, a positive slope of the best-fit regression line was used. The probabilities, P(m,ε), of finding m points within a cell of scale ε, were normalized and used to obtain the first order and higher-order mass moments: Lacunarity of a fractal set represents the ratio of the second-order moment (variance) to the first-order moment (average), of the mass probability distribution. It is a measure of structural variance within an object and an additional parameter for texture description. Higher lacunarity is related to a larger vertical spread of the regression lines corresponding to the multiple FD calculations. Denoting the mathematical expectation (average)by h�i, the fractal lacunarity was calculated as [25,26]:

Statistical analysis
Statistical analysis was performed using IBM SPSS Statistics for Windows, version 26 (IBM Corp., Armonk, N.Y., USA). The normality of the dataset was tested using Shapiro-Wilk and Kolmogorov-Smirnov tests, as well as Q-Q plots. Subsequently, image data were evaluated by analysis of variance (ANOVA) and Tukey's test for post-hoc analysis. In cases where the Levene's test showed unequal variances, we used Welch's ANOVA. Independent samples t-test was used when comparing two sets of data, as well as Mann-Whitney U for comparing data with non-normal distribution. For analysis of correlation Spearman's rank correlation test was used. The level of significance was p � 0.05.

Ethics statement
The study protocol was approved by the State Ethics Committee (Veterinary Directorate of the Ministry of Agriculture, Forestry and Water Management of Serbia decision no. 323-07-02446/2014-05/1 and 323-07-05638/2019-05). The animal experiments were conducted concordant to procedures described in the National Institutes of Health Guide for Care and Use of Laboratory Animals (Washington, DC, USA). All efforts were made to minimize suffering.

Results
We analyzed images of T. gondii cysts obtained from brain homogenates of mice infected with four different parasite strains (BGD1, BGD14, BGD26, Me49) in the same host (SW); Me49 cysts obtained at two different time points (6 and 12 weeks) in the same host (BALB/c); and Me49 cysts obtained from two different hosts (BALB/c and SW). Accordingly, the analyzed cysts differed according to the parasite strain, duration of infection (cyst age), route of infection, life stage of the parasite used for infection, and host. The total number of evaluated cysts was 31. The basic properties of the source cysts are presented in Table 1.
Parameters of interest included diameter, circularity, PD, FD, and lacunarity. A summary of all results is presented in Table 2.
The diameter of the analyzed cysts varied widely, spanning from 12.91 μm to 132.64 μm (mean 50.55 ± SD 25.37), with archived cysts being significantly larger than those obtained from the experiments (p = 0.003). Circularity was very uniform, being remarkably close to 1 for all cysts (0.993 ±0.003). Also, PD showed no statistically significant differences among cysts of different strains (p PD = 0.12), time points (p PD = 0.23), and hosts (p PD = 0.31). However, a significant negative correlation between PD and cyst diameter was observed (ρ = -0.54, p = 0.002).
Fractal analysis yielded interesting results. The FD was also very uniform among all cysts, ranging from 1.59-1.82 (mean 1.734 ± 0.063) with a variation coefficient of 0.0365 (3.6%). No difference in the FD was found between cysts of different strains (p = 0.72), time points (p = 0.51), and hosts (p = 0.57) (Fig 3A; Fig 4A and 4B). Lacunarity, on the other hand, differed among the strains, in that it was significantly lower in the BGD1 strain in comparison to the three other ones (p�0.001), as well as in BGD14 and BGD26 vs. Me49 (Tukey p = 0.044 and p = 0.033, respectively), but not between BGD14 and BGD26 (p = 0.99) (Fig 3B). Moreover, it did not differ among Me49 cysts harvested at different time points (p = 0.23), nor among those derived from different mouse strains (p = 0.56) (Fig 4C and 4D).

Discussion
Despite years of research, questions regarding the biology of T. gondii tissue cysts, such as bradyzoite organization and burden within cysts, remain unresolved. An available approach to describe and quantify the structural complexity of cells and tissues is by performing fractal analysis. The growth, branching, and rearrangement of biological entities can be considered and described as a multiplicative process. The resulting structural patterns seen under a microscope exhibit self-similarity and possess a certain degree of space-filling in comparison with other areas of the same specimen, as well as with other specimens of the same type.
In an attempt to improve the understanding of the nature of the tissue cyst structural complexity, we applied fractal analysis, as well as morphological and particle analysis, to T. gondii cysts obtained from different sources. We set out to analyze whether tissue cyst fractal geometry values are stable or tend to change depending on the parasite strain, cyst age, or mouse host genetic background. The common feature for all of the cysts was that they all originated from type-2 T. gondii strains.
Intriguingly, the majority of analyzed parameters remained uniform, i.e. circularity, PD, and FD did not differ among different parasite strains, time points, and hosts.
The analyzed cysts widely differed in diameter, with the cysts obtained from the experiment being significantly smaller, compared to historical (archived) cysts. This may be due to a bias towards larger, more interesting cysts kept in a laboratory photo archive. No other analyzed parameter (circularity, PD, FD, and lacunarity) differed between archived images and experimentally derived cysts.
According to our results, circularity was a constant parameter, with all the values very close to 1, which indicates that T. gondii cysts in the brain tend to keep an almost perfectly round shape at any stage, a result that confirms a general assumption.
FD, the value that quantifies irregularity of complex objects, can range from 1 to 2 for twodimensional images [27]. In this study, the mean FD value of all 31 analyzed cysts was 1.734 with very low variation among the cysts (SD = ± 0.063). Despite all the differences among the cysts, the FD remained stable. On the other hand, lacunarity, as a measure of heterogeneity [26], had generally moderate values (0.309 ± 0.061) but showed variation among the different parasite strains. Cysts of the BGD1 strain had significantly lower lacunarity in comparison to the three other strains, which indicates these cysts were the most homogeneous. This strain of human origin was isolated at the NRLT in 2005 and maintained in mice to date, and four out of the six analyzed images were obtained from the last, 32 nd passage. Whether the lacunarity value is an intrinsic feature of a particular strain or is related to the in vivo maintenance remains to be determined.
In the experiment performed for this study, we used the Me49 strain as a widely used reference type-2 strain. In addition to SW mice commonly used in our laboratory, we included an inbred (BALB/c) mouse strain to overcome variability associated with infection in outbred mice. Moreover, potential variability caused by poor control of the size of the infecting inoculum if infection is carried out with a specific number of cysts (containing an unknown number of bradyzoites) was overcome by using (readily countable) tachyzoites. The resulting cysts were harvested at two different time points after establishment of chronic infection to monitor potential differences among brain cysts over time i.e. potential growth of mature cysts. The results showed that, irrespective of the host genetic background and cyst age, main fractal geometry values among Me49 cysts were stable, i.e. no significant differences were found either for FD or lacunarity.
Watts et al. [7] used a different, custom-made imaging software (BradyCount 1.0) to estimate the number of bradyzoites from purified Me49 cysts based on the number of nuclei in the cysts, which were precisely visualized by the use of fluorescent dye. Calculation of PD included the watershed algorithm used for particle analysis in our study. The Watts study analyzed a large number of tissue cysts (n = 463) acquired at 5 different time points between 3 and 8 weeks post infection, and showed variations in PD, with an increase between weeks 4 and 6 but a decrease in week 8; the lower cyst occupancy at the latter point was attributed to a stabilization of the bradyzoite burden vs. the largest cyst growth. A trend of a decrease in the packing density with the increase in the cyst diameter was also observed in our study. However, we showed no difference in the PD between archived cysts (aged 6-15 months, mean ± SD = 10.71 ± 3.91) and experimentally derived cysts (aged 6-12 weeks, mean ± SD = 7.53 ± 2.29), indicating that the occupancy of the tissue cysts did not change much after 6 weeks or more post infection. This is consistent with cyst maturation; cysts have been shown to no longer grow after three months post infection [28]. The simplified methodology used here, not requiring additional manipulation of brain tissue, provided comparable results for the relationship between PD and diameter with the Watts study, indicating the potential of computational image analysis.
This study furthers the use of image analysis in describing the structural complexity of T. gondii cysts, and presents the first application of fractal analysis for this purpose. Use of a freely available software is a cost-effective mathematically quantifiable approach to advance automated image scoring for T. gondii cysts. The methodology used here is straightforward and cost-effective, requiring basic steps that are routinely applied in T. gondii bioassays. A practical application may include discriminating T. gondii cysts from similar structures including artefacts, thereby improving the specificity of bioassays [29, 30]. Moreover, as particle analysis allows for precise cell counting, it may be applied as a low cost alternative to fluorescent staining and other elaborate high-cost techniques. Since this study was limited to type-2 strains, it would be interesting to apply the same methodology to examine the structure of cysts of other T. gondii genotypes.