**3. Discussion**

In the present study a wide variety of proteins was identified and quantitated using two mass-spectrometric approaches in urine samples of patients with FSGS and MCD. No significant differences were found between the proteomic profiles of patients with MCD and a general FSGS group. However, we found that the differences in the urinary protein profiles in FSGS patients were highly dependent on the severity of the disease—thus the patients were subdivided into those with mild FSGS—normal kidney function and steroid-sensitive NS, and those characterized by severe proteinuria, impaired renal function and steroid-resistant NS. These two groups showed significant differences in the levels of 68 proteins. The group with a severe progressive FSGS showed a number of abundant kidney-derived proteins in urine.

Apolipoproteins, vitronectin, hemopexin and gelsolin reflect the process of active damage to podocytes [24–26,32,34,38–40,46–51]. Apolipoprotein A-1 (ApoA-1), a small 28 kDa high-density lipoprotein (HDL) component, is considered as a putative permeability factor [32,34,50]. In a study by Puig-Gay N. et al. an increase in ApoA-1b in the urine was noted immediately before the onset of proteinuria in patients with FSGS recurrence in the transplant [32]. This protein is absent in the urine of healthy individuals and in most patients with glomerular proteinuria caused by other glomerulopathies [34,50]. Our data also suggests a potential role of ApoA in the pathogenesis of primary FSGS. In terms of response to therapy, our data is consistent with that of Kalantari et al. who differentiated steroid-sensitive and steroid-resistant patients with a confirmed FSGS by urine proteome analysis. Among 21 proteins, ApoA-1 was one of the most significant marker between steroid-sensitive and steroid-resistant forms. An increase in ApoA-1 is associated with hyperlipidemia and low-density lipiprotein oxidation [51].

Like Kalantari et al. and other groups of researchers, we established the role of ApoA-1 and some other urinary proteins in diagnosis of FSGS [32,51,52], however, we took a different approach in the study and introduced a specific index that allows us to evaluate not only the response to steroid therapy, but also and the level of proteinuria and kidney function. This approach has been more effective in separating FSGS patients by urinary proteome to two groups. In addition, the group of patients with mild FSGS is comparable to those with MCD, while the proteomic profile of the FSGS group 2 is significantly different, and apparently, this model can be used to assess the severity of the disease.

Among the biomarkers of podocyte damage, we found an increase in the level of hemopexin, which is a glycoprotein with the highest affinity for hem [24]. Hemopexin is currently considered as one of the possible circulating factors of idiopathic FSGS. Hemopexin binds free methemoglobin, further recognized by CD91 on hepatocytes or macrophages in the spleen, liver, and bone marrow [25]. Cell culture studies have shown that hemopexin can induce the redistribution of the actin cytoskeleton in podocytes and development of proteinuria [24–26].

An increase in the urinary excretion of vitronectin in FSGS patients is also of interest, since vitronectin activates β3-integrins, molecules that ensure the fixation of podocytes to the GBM. In case of FSGS, the loss of vitronectin may be associated with the process of podocyte detachment from the glomerular basement membrane [38]. Elevated serum gelsolin aggravates the development of proteinuria and renal dysfunction by F-actin rearrangement, foot process effacement and cell movement [39,40]. These markers of podocyte dysfunction were elevated only in severe steroid-resistant FSGS.

In group 1, patients with mild FSGS, other factors of podocyte dysfunction are detected in the urine—galectin-3—binding protein. It is suggested that this protein is elevated in the kidney tissue in MCD and in the urine in lupus nephritis, and it appears to have a pro-inflammatory function [44,45].

Compared to other glomerulopathies, complement activation processes are not well understood in FSGS. In the second FSGS group we found a significant increase in urinary excretion of complement components C4b, C9, as well as factor B, I and a decrease in CD59, an inhibitor of the membrane attack complex, which indicates the possible role of complement activation associated with stronger damage. Our data is consistent with the results of a study by Huang J. et al., who showed the possibility of systemic complement activation in FSGS patients with increased levels of C3a, C5a, and C5b-9 in blood plasma and urine [43]. Activation of the alternative complement pathway in FSGS may be associated with a poor prognosis [30,31]. On the other hand, an increase in the level of some complement components in the urine may be the result of loss due to severe glomerular filter damage. For a more accurate asessement of complement system disorders in FSGS, the study these components in the blood of these patients is required.

The intensity of accumulation of ECM components is reflected by the excretion of fetuin A, β2 microglobulin, and immunoglobulin chains in the urine of patients with FSGS [37]. The same changes in the protein profiles indicating the damaging processes on tubular cells and accumulation of ECM in the interstitium can be noticed in patients with FSGS with impaired renal function resistant to steroid therapy in our study. An equally important component of disease progression and lack of response to steroid therapy is tubulointerstitial fibrosis and tubular damage in FSGS, which are reflected by an increase in the level in the urine of lumican and cystatin C. These processes are also evidenced by an increase in vitamin D-binding protein and retinol-binding protein 4 in the urine of FSGS patients with impaired renal function [18–21,37]. Simultaneously with the processes of interstitial inflammation, the mechanisms of accumulation of ECM components are activated [18–20]. Our data confirms the results of experimental studies [23,27–29,52]. A study of urinary proteins dynamic changes was conducted on a focal segmental glomerulosclerosis rat model (adriamycin-induced nephropathy) and showed that levels of Fetuin-A and alpha 1 microglobulin may be promising markers for early detection of FSGS. Thus, some proteins or their combinations can change with the disease progression [52]. Interalpha-trypsin inhibitor heavy chain H1 and H2 activate the CD44+ parietal epithelial cells, the main profibrogenic cells in the glomeruli, and thus are a powerful stimulus for glomerulosclerotic processes [27–29].

Many urine proteins are of serum origin and enter the urine through the damaged glomerular filter. However, these proteins can cause additional damage and pathology progression. For example, components of the complement system or plasminogen. The conversion of plasminogen to plasmin in urine can activate the epithelial sodium channels and cause sodium retention in the renal tubules—one of the mechanisms of renal edema [23]. Interestingly, in case of MCD and mild FSGS with a favorable prognosis, some of the proteins found in the urine seem to be of a protective nature, such as for example plasma protease C1 [35].

Perez et al. analyzed urinary peptide profiles using magnetic bead-based technology, combined with MALDI-TOF mass spectrometry, in 44 patients diagnosed with MCD and FSGS. In this work the low molecular weight fraction of urine, containing peptides of proteins were analyzed, while in our work the high-molecular protein fraction was isolated. The authors showed that FSGS patients had higher levels of uromodulin fragments and lower concentrations of fragments of A1AT [53]. In 2017 Perez V. et al. also ran a study that included 24 patients with MCD and 25 patients with FSGS and analyzed their urine proteome using two-dimensional gel electrophoresis in combination with MS to detect urinary biomarkers capable of differentiating MCD and FSGS. They showed that urine concentrations of alpha-1 antitrypsin, transferrin, histatin-3, and 39S ribosomal protein L17 were decreased in the FSGS group, and the calretinin level was increased as compared to the MCD group [54]. We did not find any significant differences between the MCD and the general FSGS groups. However, we found an increase in alpha-1 antitrypsin, and other proteins in patients with severe FSGS as compared with those with MCD and mild FSGS, who showed similar profiles of urine proteins.

Thus, we have identified a wide range of proteins that differ in patients with mild course of steroid-sensitive FSGS/MCD and FSGS patients with a progressive steroidresistant NS. Proteins excreted in the urine reflect damage to several parts of the nephron podocytes, tubulo-interstitium, accumulation of ECM, complement activation. Patients

with MCD and FSGS with steroid-sensitivity and a favorable course of progress showed similar urine protein profiles, while severe progressive FSGS with steroid-resistant nephrotic syndrome were characterized by the early activation of the complement system and profibrogenic mechanisms and accumulation of extracellular matrix components at the onset of the disease.


**2.1 Governing equations for transport of scalars: Unclosed statistical set**

The turbulent transfer equations for a scalar *F* are usually expressed as

equation furnished an upper limit for the normalized RMS value, which is not reached

The random square waves were also used by Schulz et al. (1991) to quantify the so called "intensity of segregation" in the superficial boundary layer formed during mass transport, for which the explanations of segregation scales found in Brodkey (1967) were used. The time constant of the intensity of segregation, as defined in the classical studies of Corrsin (1957, 1964), was used to correlate the mass transfer coefficient across the water surface with more usual parameters, like the Schmidt number and the energy dissipation rate. Random square waves were also applied by Janzen (2006), who used the techniques of Particle Image Velocimetry (PIV) and Laser Induced Fluorescence (LIF) to study the mass transfer at the air-water interface, and compared his measurements with the predictions of Schulz & Schulz (1991) employing *ad hoc* concentration profiles. Further, Schulz & Janzen (2009) confirmed the upper limit for the normalized RMS of the concentration fluctuations by taking into account the effect of diffusion, also evaluating the thickness of diffusive layers and the role of diffusive and turbulent transports in boundary layers. A more detailed theoretical relationship for the RMS of the concentration fluctuation showed that several

Intending to present the methodology in a more organized manner, Schulz et al. (2011a) showed a way to "model" the records of velocity and mass concentration (that is, to represent them in an *a priori* simplified form) for a problem of mass transport at gas-liquid interfaces. The fluctuations of these variables were expressed through the so called "partition, reduction, and superposition functions", which were defined to simplify the oscillating records. As a consequence, a finite number of basic parameters was used to express all the statistical quantities of the equations of the problem in question. The extension of this approximation to different Transport Phenomena equations is demonstrated in the present study, in which the mentioned statistical functions are derived for general scalar transport (called here "scalar-velocity interactions"). A first application for velocity fields is also shown (called here "velocity-velocity interactions"). A useful consequence of this methodology is that it allows to "close" the turbulence equations, because the number of equations is bounded by the number of basic parameters used. In this chapter we show 1) the *a priori* modeling (simplified representation) of the records of turbulent variables, presenting the basic definitions used in the random square wave approximation (following Schulz et al., 2011a); 2) the generation of the usual statistical quantities considering the random square wave approximation (scalar-velocity interactions); 3) the application of the methodology to a one-dimensional scalar transport problem, generating a closed set of equations easy to be solved with simple numerical resources; and 4) the extension of the study of Schulz & Johannes (2009) to velocity fields (velocity-velocity

Because the method considers primarily the oscillatory records itself (*a priori* analysis), and not phenomenological aspects related to physical peculiarities (*a posteriori* analysis, like the definition of a turbulent viscosity and the use of turbulent kinetic energy and its dissipation

rate), it is applicable to any phenomenon with oscillatory characteristics.

**2.1 Governing equations for transport of scalars: Unclosed statistical set**  The turbulent transfer equations for a scalar *F* are usually expressed as

different statistical profiles of turbulent mass transfer may be interrelated.

when diffusion is taken into account.

interactions).

**2. Scalar-velocity interactions**

$$\frac{\partial \overline{\mathbf{F}}}{\partial \mathbf{\hat{t}}} + \overline{\mathbf{V}}\_{\mathbf{i}} \frac{\partial \overline{\mathbf{F}}}{\partial \mathbf{x}\_{\mathbf{i}}} = \frac{\partial}{\partial \mathbf{x}\_{\mathbf{i}}} \left( \mathbf{D}\_{\mathbf{F}} \frac{\partial \overline{\mathbf{F}}}{\partial \mathbf{x}\_{\mathbf{i}}} - \overline{\mathbf{v}\_{\mathbf{i}} \mathbf{f}} \right) + \overline{\mathbf{g}} \ \mathbf{y} \ \mathbf{i} = \mathbf{1} \ \mathbf{2} \ \mathbf{3} \tag{1}$$

where *F* and *f* are the mean scalar function and the scalar fluctuation, respectively. *Vi* (*i* = 1, 2, 3) are mean velocities and *vi* are velocity fluctuations, *t* is the time, *xi* are the Cartesian coordinates, *g* represents the scalar sources and sinks and *DF* is the diffusivity coefficient of *F*. For one-dimensional transfer, without mean movements and generation/consumption of *F*, equation (1) with *x3*=*z* and *v3*=߱ is simplified to

$$\frac{\partial \overline{\mathbf{F}}}{\partial \mathbf{f}} = \frac{\partial}{\partial \mathbf{z}} \left( \mathbf{D}\_{\mathbf{F}} \frac{\partial \overline{\mathbf{F}}}{\partial \mathbf{z}} - \overline{\mathbf{a} \, \mathbf{f}} \right) \tag{2}$$

As can be seen, a second variable, given by the mean product ωf , is added to the equation of *F* , so that a second equation involving ωf and *F* is needed to obtain solutions for both variables. Additional statistical equations may be generated averaging the product between equation (1) and the instantaneous fluctuations elevated to some power ( *f* ). As any new equation adds new unknown statistical products to the problem, the resulting system is never closed, so that no complete solution is obtained following strictly statistical procedures (closure problem). Studies on turbulence consider a low number of statistical equations (involving only the first statistical moments), together with additional equations based on *ad hoc* models that close the systems. This procedure seems to be the most natural choice, because having already obtained equation (2), it remains to model the new parcel *ω f a posteriori* (that is, introducing hypotheses and definitions to solve it). An example is the combined use of the Boussinesq hypothesis (in which the turbulent viscosity/diffusivity is defined) with the Komogoroff reasoning about the relevance of the turbulent kinetic energy and its dissipation rate. The ߢെߝ model for statistical turbulence is then obtained, for which two new statistical equations are generated, one of them for *k* and the other for ߝ. Of course, new unknown parameters appear, but also additional *ad hoc* considerations are made, relating them to already defined variables.

In the present chapter, as done by Schulz et al. (2011a), we do not limit the number of statistical equations based on *a posteriori* definitions for *ω f* . Convenient *a priori* definitions are used on the oscillatory records, obtaining transformed equations for equation (1) and additional equations. The central moments of the scalar fluctuations, *f F F* , ߠ =1, 2, 3,… are considered here. For example, the one-dimensional equations for ߠ=2, 3 and 4, are given by

$$\frac{1}{2}\frac{\partial \overline{f^2}}{\partial t} + \overline{f}\overline{o}\frac{\partial \overline{F}}{\partial z} + \frac{1}{2}\frac{\partial \overline{f^2}}{\partial z} = D\_F \left( \overline{f \frac{\partial^2 f}{\partial z^2}} \right) \tag{3a}$$

$$\frac{1}{3}\frac{\overline{\partial}\overline{f^3}}{\overline{\partial}t} + \overline{f^2}\frac{\partial\overline{F}}{\overline{\partial}t} + \overline{f^2}o\frac{\partial\overline{F}}{\overline{\partial}z} + \frac{1}{3}\frac{\overline{\partial}o\overline{f^3}}{\overline{\partial}z} = D\_\mathbb{F}\left(\overline{f^2}\frac{\partial^2\overline{F}}{\overline{\partial}z^2} + \overline{f^2}\frac{\overline{\partial^2f}}{\overline{\partial}z^2}\right) \tag{3b}$$

One Dimensional Turbulent Transfer

**field -** *n*

depend on *z*.

for example, the central moments *f F F*

oscillations justifies these corrective parcels.

This definition also implies that

Using Random Square Waves – Scalar/Velocity and Velocity/Velocity Interactions 7

2 1 *tt t* is the time interval for the average operation. Equation (4) generates a mean value *F z*( ) for 1 2 *z zz* and 1 2 *t tt* . Any statistical quantity present in equations 3, like,

The method described in the next sections allows to obtain the relevant statistical quantities

**2.3 Bimodal square wave: Mean values using a time-partition function for the scalar**

The basic assumptions made to "model" the original oscillatory records may be followed considering Figure 2. In this sense, figure 2a is a sketch of the original record of the scalar variable *F* at a position 1 2 *z zz* , as shown in the gray vertical plane of Figure 1. The objective of this analysis is to obtain an equation for the mean function *F z*( ) for 1 2 *t tt* , which is also shown in figure 2a. The values of the scalar variable during the turbulent transfer are affected by both the advective turbulent movements and diffusion. Discarding diffusion, the value of *F* would ideally alternate between the limits *Fp* and *Fn* (the bimodal square wave), as shown in Figure 2b (the fluid particles would transport only the two mentioned *F* values). This condition was assumed as a first simplification, but maintaining the correct mean, in which *F z*( ) is unchanged. It is known that diffusion induces fluxes governed by *F* differences between two regions of the fluid (like the Fourier law for heat transfer and the Fick law for mass transfer). These fluxes may significantly lower the amplitude of the oscillations in small patches of fluid, and are taken into account using *Fp*-*P* and *Fn*+*N* for the two new limiting *F* values, as shown in Figure 2c. The parcels *P* and *N*

In other words, the amplitude of the square oscillations is "adjusted" (modeled), in order to approximate it to the mean amplitude of the original record. As can be seen, the aim of the method is not only to evaluate *F* adequately, but also the lower order statistical quantities that depend on the fluctuations, which are relevant to close the statistical equations. The parcels *P* and *N* were introduced based on diffusion effects, but any cause that inhibits

The first statistical parameter is represented by *n*, and is defined as the fraction of the time for which the system is at each of the two *F* values (equations 5 and 6), being thus named as

> ( ) of the observation *p t at F P*

> > of the observation *n t at F N*

(5)

(6)

"partition function". This function n depends on *z* and is mathematically defined as

( ) 1

*F* remains the same in figures 2a, b and c. The constancy between figures 2b and c is obtained using mass conservation, implying that *P* and *N* are related through equation (7):

*t*

*n*

*n*

*t*

simplify notation, both coordinates (*z*, *t*) are dropped off in the rest of the text.

of the governing equations, like the mean function *F* , using simplified records of *F*.

, is defined according to equation (4). To

$$\frac{1}{4}\frac{\partial\overline{f^4}}{\partial t} + \overline{f^3}\frac{\partial\overline{F}}{\partial t} + \overline{f^3}o\frac{\partial\overline{F}}{\partial z} + \frac{1}{4}\frac{\partial\overline{o\,f^4}}{\partial z} = D\_{\overline{F}}\left(\overline{f^3}\frac{\partial^2\overline{F}}{\partial z^2} + \overline{f^3}\frac{\partial^2\overline{f}}{\partial z^2}\right) \tag{3c}$$

In this example, equation (3a) involves *F* and *f* of equation (2), but adds three new unknowns. The first four equations (2) and (3 a, b, c) already involve eleven different statistical quantities: *F* , 2 *f* , 3 *f* , 4 *f* , *f* , 2 *f* , 3 *f* , 4 *f* , 2 2 *f f z* , 2 2 2 *f f z* , and

2 3 2 *f f z* , and the "closure" is not possible. The general equation for central moments, for

any ߠ, is given by [20]

$$\frac{1}{\partial} \frac{\overline{\partial \, f^{\theta}}}{\partial t} + \overline{f^{\theta - 1}} \frac{\partial \overline{F}}{\partial t} + \overline{f^{\theta - 1}} \alpha \frac{\partial \overline{F}}{\partial z} + \frac{1}{\theta} \frac{\partial \overline{\alpha \, f^{\theta}}}{\partial z} = D\_{\overline{F}} \left( \overline{f^{\theta - 1}} \frac{\partial^{2} \overline{F}}{\partial z^{2}} + \overline{f^{\theta - 1}} \frac{\overline{\partial^{2} f}}{\partial z^{2}} \right) \tag{3d}$$

(using ߠ=1 reproduces equation (2)).

As mentioned, the method models the records of the oscillatory variables, using random square waves. The number of equations is limited by the number of the basic parameters defined "*a priori*".

### **2.2 "Modeling" the records of the oscillatory variables**

As mentioned in the introduction, the term "modeling" is used here as "representing in a simplified way". Following Schulz et al. (2011a), consider the function *F*(*z*, *t*) shown in Figure 1. It represents a region of a turbulent fluid in which the scalar quantity *F* oscillates between two functions *Fp* (*p*=previous) and *Fn* (*n*=next) in the interval *z*1<*z*<*z*2. Turbulence is assumed stationary.

Fig. 1a. A two-dimensional random scalar field F oscillating between the boundary functions *Fp*(t) and *Fn*(t).

Fig. 1b. Sketch of the region shown in figure (1a). Turbulence is stationary.  The first four equations (2) and (3 a, b, c) already involve eleven different

 , 2 *f* , 3 *f* , 4 *f* , 2

, and the "closure" is not possible. The general equation for central moments, for

1 1 1 1

*f F F f F f f f Df f t t zz z z*

As mentioned, the method models the records of the oscillatory variables, using random square waves. The number of equations is limited by the number of the basic parameters

As mentioned in the introduction, the term "modeling" is used here as "representing in a simplified way". Following Schulz et al. (2011a), consider the function *F*(*z*, *t*) shown in Figure 1. It represents a region of a turbulent fluid in which the scalar quantity *F* oscillates between two functions *Fp* (*p*=previous) and *Fn* (*n*=next) in the interval *z*1<*z*<*z*2. Turbulence is

*F*

1 1

1 1

**2.2 "Modeling" the records of the oscillatory variables**

In this example, equation (3a) involves *F* and *f*

statistical quantities: *F* , 2 *f* , 3 *f* , 4 *f* , *f*

(using ߠ=1 reproduces equation (2)).

Fig. 1a. A two-dimensional random scalar field F oscillating between the boundary

functions *Fp*(t) and *Fn*(t).

2 3 2 *f f z*

any ߠ, is given by [20]

defined "*a priori*".

assumed stationary.

4 4 *F*

2 2

2 2

Fig. 1b. Sketch of the region shown in figure (1a). Turbulence is stationary. Adapted from

1 2

(4)

Schulz et al. (2011a)

The time average of *F(z, t)* for 1 2 *z zz* , indicated by *Fzt* (,) is defined as usual

*F z t F z t dt for z z z t*

2

1 (, ) (,) *t*

1

*t*

2 2

of equation (2), but adds three new

2 *f f z*

 , 2 2 2 *f f z*

(3c)

, and

(3d)

2 1 *tt t* is the time interval for the average operation. Equation (4) generates a mean value *F z*( ) for 1 2 *z zz* and 1 2 *t tt* . Any statistical quantity present in equations 3, like, for example, the central moments *f F F* , is defined according to equation (4). To

simplify notation, both coordinates (*z*, *t*) are dropped off in the rest of the text.

The method described in the next sections allows to obtain the relevant statistical quantities of the governing equations, like the mean function *F* , using simplified records of *F*.


*Article* **A Shock-in-Jet Synchrotron Mirror Model for Blazars**

**Markus Böttcher**

Centre for Space Research, North-West University, Potchefstroom 2520, South Africa<EMAIL_ADDRESS>
**Abstract:** Reinhard Schlickeiser has made groundbreaking contributions to various aspects of blazar physics, including diffusive shock acceleration, the theory of synchrotron radiation, the production of gamma-rays through Compton scattering in various astrophysical sources, etc. This paper, describing the development of a self-consistent shock-in-jet model for blazars with a synchrotron mirror feature, is therefore an appropriate contribution to a Special Issue in honor of Reinhard Schlickeiser's 70th birthday. The model is based on our previous development of a self-consistent shock-in-jet model with relativistic thermal and non-thermal particle distributions evaluated via Monte-Carlo simulations of diffusive shock acceleration, and time-dependent radiative transport. This model has been very successful in modeling spectral variability patterns of several blazars, but has difficulties describing orphan flares, i.e., high-energy flares without a significant counterpart in the low-frequency (synchrotron) radiation component. As a solution, this paper investigates the possibility of a synchrotron mirror component within the shock-in-jet model. It is demonstrated that orphan flares result naturally in this scenario. The model's applicability to a recently observed orphan gamma-ray flare in the blazar 3C279 is discussed and it is found that only orphan flares with mild ( a factor of 2–3) enhancements of the Compton dominance can be reproduced in a synchrotronmirror scenario, if no additional parameter changes are invoked.

**Keywords:** active galaxies; blazars; diffusive shock acceleration; radiative transport; gamma-rays

## **1. Introduction**

Blazars are a class of jet-dominated active galactic nuclei. As most convincingly argued by Reinhard Schlickeiser (RS) in 1996 [1], their broad-band non-thermal emission, ranging from radio to gamma-rays, must be strongly Doppler boosted due to relativistic motion of an emission region along the jet, oriented close to our line of sight. The spectral energy distributions (SEDs) of blazars are dominated by two broad, non-thermal radiation components. The low-frequency component, from radio to optical/UV/X-ray frequencies, is generally attributed to synchrotron radiation by relativistic electrons. Most notably, Crusius and Schlickeiser [2,3] have evaluated the angle-averaged synchrotron emission from isotropically distributed electrons in random magnetic fields, including plasma effects, which are now frequently used as the standard expressions for the low-frequency emission from blazars. However, note also an alternative suggestion by RS in 2003 [4] that the lowfrequency emission may be produced as electrostatic bremsstrahlung, i.e., the scattering of electrostatic Langmuir waves excited by two-stream instabilities, as expected in the jet-inter-stellar-medium interaction scenario of Schlickeiser et al. (2002) [5].

Motivated by early *γ*-ray observations by the SAS-2 and COS-B satellites, already in 1979–1980, RS had considered inverse-Compton scattering as the dominant mechanism to produce high-energy *γ*-rays in astrophysical sources, pointing out the importance of Klein-Nishina effects in the calculation of *γ*-ray spectra [6–8]. Also in leptonic models for blazars, inverse-Compton scattering by relativistic electrons in the jet is considered the dominant high-energy emission mechanism. Target photons for Compton scattering can be the co-spatially produced synchrotron (or electrostatic bremsstrahlung) radiation, in which case it is termed synchrotron self-Compton (SSC) emission (e.g., [9,10]). The first suggestion of target photon fields from outside the jet involved RS in two seminal papers suggesting

**Citation:** Böttcher, M. A Shock-in-Jet Synchrotron Mirror Model for Blazars. *Physics* **2021**, *3*, 1112–1122. https://doi.org/10.3390/ physics3040070

Received: 16 September 2021 Accepted: 5 November 2021 Published: 22 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the photon field of the accretion disk as the dominant target photon field [11,12]. Alternative sources of external target photons may be the broad-line region (BLR) (e.g., [13]), a dusty, infra-red emitting torus (e.g., [14]), or other regions of the jet (e.g., [15,16]). The relativistic motion of the high-energy emission region in a blazar jet through these generally anisotropic external radiation fields leads to complicated transformation properties from the active galactic nucleus (AGN) rest frame into the emission-region frame, which were studied in detail by Dermer and Schlickeiser in 2002 [17]. Which of these potential radiation fields might dominate, depends critically on the location of the emission region, which can be constrained by the absence of obvious signatures of *γγ* absorption of high-energy and very-high-energy *γ*-rays by the nuclear radiation fields of the central AGN, with one of the first detailed discussions of such constraints published by Dermer and Schlickeiser in 1994 [18].

The generation of the non-thermal broadband emission from blazars requires the efficient acceleration of electrons to ultra-relativistic energies. One of the plausible mechanisms of particle acceleration acting in the relativistic jets of blazars is diffusive shock acceleration (DSA), which was studied in the context of a general derivation of the kinetic equation of test particles in turbulent plasmas by RS in two seminal papers in 1989 [19,20] for nonrelativistic shock speeds, while particle acceleration by magnetic turbulence, specifically in relativistic jets was studied by Schlickeiser and Dermer in 2000 [21]. Particle acceleration at relativistic shocks has been considered by several authors, using both analytical methods (e.g., [22–24]) and Monte-Carlo techniques (e.g., [25–29]). The simulations by Niemiec and Ostrowski [28] and Summerlin and Baring [29] indicate that diffusive shock acceleration at oblique, mildly relativistic shocks is able to produce relativistic, non-thermal particle spectra with a wide range of spectral indices, including as hard as *n*(*p*) ∝ *p*−1, where *p* is the particle's momentum.

In two recent papers [30,31], we had coupled Monte-Carlo simulations of diffusive shock acceleration (DSA), using the code of Summerlin and Baring [29], with timedependent radiation transfer, based on radiation modules originally developed by Böttcher, Mause and Schlickeiser in 1997 [32] and further developed as detailed in [33,34]. In those studies, we found that the particles' mean free path for pitch-angle scattering, *λ*pas, which mediates the first-order Fermi process in DSA, must have a strong dependence on particle momentum, with an index *α* > 1 for a parameterization of *λ*pas(*p*) ∝ *pα*. This likely indicates a decaying level of magneto-hydrodynamic turbulence with increasing distance from the shock front. Higher-energy particles, with their larger gyro radii, then probe more distant regions from the shock front, experiencing less efficient pitch-angle scattering. Time-dependent simulations of DSA plus radiation transfer were used to fit the multiwavelength variability of the blazars 3C279 and Mrk 501 in [31] and the X-ray variability of 1ES 1959 + 650 in [35]. Multi-wavelength flares with approximately equal flare amplitude in the low-frequency (synchrotron) and high-frequency (Compton) components of the SED were naturally produced by an increase of the power injected into shock-accelerated particles, without the need for significant changes of the plasma parameters determining *λ*pas(*p*).

However, an orphan *γ*-ray flare on 20 December 2013, with no significant counterpart in the synchrotron emission component, reported as Flare B in [36], presented a severe challenge to this as well as any other single-zone emission model for blazars. A fit to the observed *γ*-ray flare was possible with a significant hardening of the DSA-generated particle spectrum as the result of a reduction of the pitch-angle-scattering mean-free path, both in overall normalization *λ*pas(0) and index *α*. However, keeping the optical (synchrotron) flux approximately constant, as observed, required a reduction of the magnetic field by a factor of 8.7, followed by a gradual recovery to the quiescent-state value with a fine-tuned time dependence. While the authors argue that such magnetic-field reductions and subsequent gradual recoveries after the passage of a shock have indeed been observed in interplanetary shocks (e.g., [37]), it is worth exploring alternative ways to explain orphan *γ*-ray flares in blazars within the framework of the shock-in-jet model developed in [30,31].

One plausible way of producing orphan *γ*-ray flares in the framework of a leptonic single-zone blazar model is the temporary enhancement of an external radiation field that serves as target for inverse-Compton scattering. This is the basis of a class of models termed *synchrotron mirror models*, where the synchrotron radiation of the high-energy emission region traveling along the jet, is reflected by a cloud to re-enter the emission region at a later time. Such models were first considered by Ghisellini and Madau [38], however without proper consideration of light-travel time effects, and by Böttcher and Dermer [39] and Bednarek [40], properly treating light-travel time effects, but considering primarily the time-dependence of the target-photon energy density without detailed calculations of the emerging *γ*-ray spectra. The synchrotron mirror model was more recently re-visited by Vittorini et al. [41], with a fully time-dependent leptonic synchrotron mirror model applied to the spectral variability of 3C454.3 in 2010 November, and Tavani et al. [42], considering also moving mirrors and applying the model to the light curve of the same flare B of 3C279 considered by [31]. Note a similar model termed the "ring of fire" model by MacDonald et al. [43,44], where the emission region passes a static synchrotron-emitting region of an outer sheath of the jet (the "ring of fire"), which produces very similar variability features as the synchrotron mirror model.

In the present paper, the time-dependent shock-in-jet model of Böttcher and Baring [31] is extended to include self-consistently a synchrotron-mirror component. Section 2 describes the additions to the model. Section 3 presents the resulting spectral variability features from an attempt to apply this model to the orphan *γ*-ray flare B of 3C279. Section 4 summarizes and discusses the results.

#### **2. Model Description**

The model developed here is a further development of the time-dependent shock-injet model of Böttcher and Baring [31]. In addition to the radiation components already included in [31] , we now introduce synchrotron emission reflected by a spherical cloud of radius *R*cl at a distance *z*cl from the central engine, assumed for simplicity to be located close to the path of the jet, however, not hydrodynamically interacting with it, as considered, e.g., by the jet-star interaction model [45,46] or the cloud ablation model [47,48]. A mildly relativistic, oblique shock is propagating along the jet, thus accelerating particles in the local environment of the shock which constitutes our moving emission region of radius *Rb*. The emission region is starting out at time *te* = 0 (in the AGN rest frame) at a height *z*0 above the black-hole—accretion-disk system powering the jet, and is propagating with a bulk Lorentz factor Γ, corresponding to a speed of *β*Γ*c*. Thus, at any given time *te*, the emission region is located at *ze* = *z*0 + *β*Γ *c te*.

Synchrotron radiation emitted by the emission region at *ze* is reflected back by the cloud to re-enter the emission region at a distance *zr* from the central engine, given by

$$z\_{\tau} = \frac{2\beta\_{\Gamma}z\_{\mathrm{cl}} + z\_{\mathrm{e}}(1 - \beta\_{\Gamma})}{1 + \beta\_{\Gamma}},\tag{1}$$

at a time (in the AGN rest frame) *tr* given by

$$t\_r = t\_\mathfrak{e} + 2\frac{z\_{\mathfrak{e}1} - z\_\mathfrak{e}}{(1 + \beta\_\Gamma)\mathfrak{e}}.\tag{2}$$

Equation (2) may be inverted to find the time at which reflected synchrotron radiation received at time *tr* has been emitted:

$$t\_{\varepsilon} = \frac{1 + \beta\_{\Gamma}}{1 - \beta\_{\Gamma}} t\_r - \frac{2}{1 - \beta\_{\Gamma}} \frac{z\_{\rm cl} - z\_0}{c} \tag{3}$$

implying that reflected synchrotron emission will be received starting at a time *t*0 (corresponding to *te* = 0) given by

$$t\_0 = \frac{2}{1 + \beta\_\Gamma} \frac{z\_{\rm cl} - z\_0}{c}. \tag{4}$$

Reflected synchrotron radiation will be received by the emission region until it passes the cloud at *t*pass ≈ (*z*cl − *z*0)/(*β*Γ*c*).

The code writes out the observed synchrotron emission spectra, *νF*sy *ν* (*te*) for every time step (with times in the observer's frame) as the shock propagates along the jet. Therefore, at any time *t*AGN > *t*0, one can use Equation (3) to find the time (in the AGN frame) at which synchrotron radiation reflected back into the emission region, has been emitted. The synchrotorn flux irradiating the cloud, *νF*cl *ν* , is then found as

$$\nu F\_{\nu}^{\rm cl} = \nu F\_{\nu}^{\rm sy}(t\_{\varepsilon}) \frac{d\_L^2}{(z\_{\rm cl} - z\_{\rm c})^2} \tag{5}$$

where *dL* is the luminosity distance to the source.

Assuming, for simplicity, that the cloud re-radiates a fraction *τ*cl of the impinging synchroton radiation isotropically, it will emit a spectral luminosity of *νL*cl nu = *πR*2 cl *τ*cl *νF*cl *ν* . The emission region will thus receive a flux of Reflected Synchrotron (RS—happy coincidence) radiation, in the comoving frame, of

$$\nu' F\_{\nu'}^{\rm RS}(t\_r) \approx \frac{R\_{\rm cl}^2 \,\mathrm{r\_{cl}} \,\nu F\_{\nu}^{\rm sy}(t\_r) \Gamma^4 \,d\_L^2}{4 \,(z\_{\rm cl} - z\_{\rm c})^2 \,(z\_{\rm cl} - z\_r)^2} \tag{6}$$

where *ν* ≈ Γ*ν* is the photon frequency in the co-moving frame.

The code evaluates a time-evolving reflected-synchrotron photon field in the emission region, *n* RS( , *t r*), where = *hν* /(*mec*2) is the dimensionless photon energy in the emission-region frame by the interplay of RS emission entering the emission region at a rate *dn* RS,inj( , *t r*)/*dt r* = *π R*2 *bν F*RS *ν* (*tr*)/(*Vb* 2 *mec*2), where *Vb* is the volume of the emission region, and escape on an escape time scale *t* esc = 3 *Rb*/(4 *c*), over a simulation time step Δ*t* as

$$
\Delta n\_{\rm RS}^{\prime}(\boldsymbol{\epsilon}^{\prime}, \boldsymbol{t}\_{r}^{\prime}) = \left(\frac{\pi R\_{\rm b}^{2}}{V\_{\rm b}} \frac{\boldsymbol{\nu}^{\prime} F\_{\boldsymbol{\nu}^{\prime}}^{\rm RS}(t\_{r})}{\boldsymbol{\epsilon}^{\prime 2} m\_{\boldsymbol{\epsilon}} \boldsymbol{\varepsilon}^{2}} - \frac{n\_{\rm RS}^{\prime}(\boldsymbol{\epsilon}^{\prime}, \boldsymbol{t}\_{r}^{\prime})}{t\_{\rm csc}^{\prime}}\right) \Delta t^{\prime}. \tag{7}
$$

In the above expressions, *h* is the Planck constant, *me* the electron mass, and *t r* = *tr*/Γ.

The time-dependent RS photon field resulting from Equation (7) acts target for inverse-Compton scattering to produce the synchrotron-mirror Compton emission, and synchrotronmirror Compton cooling is included self-consistently. For the evaluation of the synchrotronmirror Compton emission, it is assumed, for simplicity, that the target photons enter the jet directly from the front, and the head-on approximation (e.g., [49]) for the Compton cross section is used. Hence, photons scattered along the viewing direction, making an angle *θ* obs with respect to the jet axis in the co-moving frame of the emission region, with *μ* obs ≡ cos *θ* obs, have been scattered by a scattering angle *μ* = −*μ* obs. The rate density at which Reflected Synchrotron Compton (RSC) emission is produced, is calculated as

$$\dot{n}'\_{\rm RSC}(\boldsymbol{\epsilon}'\_{\rm s}, \boldsymbol{\Omega}'\_{\rm s}) = c \int\_0^\infty d\boldsymbol{\epsilon}' \, n'\_{\rm RS}(\boldsymbol{\epsilon}') \int\_1^\infty d\boldsymbol{\gamma} \, \left(1 + \beta \, \mu'\_{\rm obs}\right) n\_\boldsymbol{\epsilon}(\boldsymbol{\gamma}) \, \frac{d\sigma\_\mathbb{C}}{d\boldsymbol{\epsilon}'\_{\rm s}}(\boldsymbol{\epsilon}', \boldsymbol{\epsilon}'\_{\rm s}, \mu') \tag{8}$$

where

$$\frac{d\sigma\_{\mathbb{C}}}{d\epsilon\_{\mathbb{s}}^{\prime}} = \frac{\pi}{\gamma} \frac{r\_{\varepsilon}^{2}}{\epsilon\_{\mathbb{s}}} \left\{ y + \frac{1}{y} - \frac{2\,\epsilon\_{\mathbb{s}}^{\prime}}{\gamma\,\epsilon\_{\mathbb{c}}\,y} + \left(\frac{\epsilon\_{\mathbb{s}}^{\prime}}{\gamma\,\epsilon\_{\mathbb{c}}\,y}\right)^{2} \right\} \tag{9}$$

with *e* = *γ* (1 + *β μ* obs), *y* = 1 − *s*/*γ* and *β* = 1 − 1/*γ*2 [49].

#### **3. Results: Spectral Variability Features**

As detailed in the introduction, the study of the synchrotron mirror model developed here was motivated by the difficulties in modeling the orphan *γ*-ray flare B of 3C279 in December 2013 reported by Hayashida et al. [36]. We therefore start with the quiescent-state parameters of the shock-in-jet model for 3C279 used in [31]. The emission region is set to start out at *z*0 = 0.1 pc, and the cloud acting as the mirror is assumed to be located at *z*cl = 1 pc. The bulk Lorentz factor is Γ = 15, as used in [31]. The cloud radius is assumed to be *R*cl = 3 × 1017 cm and its reflective fraction is *τ*cl = 0.1. The complete list of model parameters can be found in Table 1.

**Table 1.** Relevant model parameters for the case study motivated by the December 2013 flare of 3C279.

The resulting sequence of snap-shot SEDs (starting right before the onset of the synchrotron-mirror Compton emission) is illustrated in Figure 1. It is clear that the model does produce a significant orphan *γ*-ray flare, accompanied by a slight reduction of the synchrotron emission due to the increased Compton cooling of relativistic electrons. The latter is consistent with the observed evolution of the SED. However, the amplitude of the orphan *γ*-ray flare amounts only to an enhancement of the *γ*-ray flux by a factor of ∼2, in contrast to the observed dramatic increase by a factor of 10. Even increasing the synchrotron-mirror efficiency (e.g., by increasing *τ*cl or *R*cl; see Equation (6)), does not increase the *γ*-ray flare amplitude substantially. The reason for this is that the *γ*-ray flux is limited by the available power injected into shock-accelerated electrons, which are already in the fast-cooling regime, thus radiating very efficiently. Any further enhancement of the reflected-synchrotron energy density will only suppress the synchrotron emission further, but not lead to a significant increase of the *γ*-ray flare amplitude. We therefore conclude that a pure shock-in-jet synchrotron mirror scenario is not able to produce the observed large-amplitude orphan *γ*-ray flare in 3C279 in December 2013. In order to achieve this, additional power would need to be injected into shock-accelerated electrons, leaving us with the same difficulties encountered in [31], i.e., requiring a fine-tuned reduction and gradual recovery of the magnetic field.

Nevertheless, in spite of its inapplicability to this particular orphan flare, it is worthwhile considering this simulation for a generic study of the expected spectral variability patterns in the shock-in-jet synchrotron mirror model. The multi-wavelength light curves at 5 representative frequencies (high-frequency radio, optical, X-rays, high-energy [HE, 200 MeV], and very-high-energy [VHE, 200 GeV] *γ*-rays) are shown in Figure 2. All light curves in the Compton SED component (X-rays to VHE *γ*-rays) show a flare due to the synchrotron-mirror Compton emission. Note that the VHE *γ*-ray light curve had to be scaled up by a factor of 1010 to be visible on this plot. Thus, the apparently large VHE flare is actually at undetectably low flux levels for the parameters chosen here. In contrast, the 230 GHz radio and optical light curves show a dip due to increased radiative cooling during the synchrotron mirror action. The radio dip is significantly delayed compared to the optical due to the longer cooling time scales of electrons emitting in the radio band.

**Figure 1.** Spectral energy distributions (SEDs) of 3C279 in 2013–2014, from [36], along with snap-shot model SEDs from the shock-in-jet synchrotron-mirror model. The dashed vertical lines indicate the frequencies at which light curves and hardness-intensity relations were extracted. The legend follows the nomenclature of different periods from Hayashida et al. (2015) [36].

**Figure 2.** Model light curves in various frequency/energy bands resulting from the synchrotron mirror simulation illustrated in Figure 1 at the 5 representative frequencies/energies marked by the vertical dashed lines. Note that the very-high-energy (VHE, 200 GeV) *γ*-ray flux is scaled up by a factor of 1010 in order to be visible on the plot.

Cross-correlation functions between the various light curves from Figure 2 are shown in Figure 3. As expected from inspection of the light curves, significant positive correlations between X-rays and the 2 *γ*-ray bands with only small time lags (*γ*-rays leading X-rays by a few hours) and between the radio and optical band, with optical leading the radio by ∼15 h, are seen. The synchrotron (radio and optical) light curves are anti-correlated with the Compton (X-rays and *γ*-rays) ones, again with a significant lag of the radio emission by ∼15 h.

**Figure 3.** Cross-correlation functions between the model light curves in various energy/frequency bands.

Figure 4 shows the hardness-intensity diagrams for the 5 selected frequencies/energies, i.e., the evolution of the local spectral index (*a*, defined by *Fν* ∝ *ν*−*a*) vs. differential flux. Generally, all bands, except the optical, exhibit the frequently observed harder-whenbrighter trend. Only the radio and X-ray bands show very moderate spectral hysteresis. The dip in the optical R-band) light curve is accompanied by a very slight redder-whenbrighter trend, likely as the consequence of the increased radiative cooling during the synchrotron-mirror action.

**Figure 4.** Model hardness-intensity diagrams at the selected frequencies/energies. The spectral index *a* is defined by *Fν* ∝ *ν*−*a*, so that a smaller value indicates a harder spectrum. The VHE band has been omitted here due to its unobservably low flux level and very steep local spectral index. Arrow indicate the evolution in time.

#### **4. Summary, Discussion, and Conclusions**

In this paper, the leptonic shock-in-jet blazar model of [31] is extended with the addition of a self-consistent synchrotron mirror component. This was motivated by the difficulty in modeling orphan *γ*-ray flares with such an effectively single-zone model. A particularly high-amplitude (factor of ∼10) orphan *γ*-ray flare of the blazar 3C279 from December 2013 was chosen as a case study. However, the attempt to model this flare with the shock-in-jet synchrotron mirror model developed here, failed because the maximum *γ*-ray flux was limited by the (fixed) amount of power injected into shock-accelerated electrons, allowing for orphan flares with amplitudes of at most ∼2–3. Higher-amplitude flares would require an enhanced energy injection into relativistic electrons, in addition to more efficient pitch-angle scattering, leading to a harder electron spectrum. However, this would cause the same difficulties of having to decrease the magnetic field, followed by a fine-tuned recovery to its quiescent state, as were encountered in [31].

More successful model representations of this particular flare of 3C279 were presented by several authors. Hayashida et al. [36] use the model of Nalewajko et al. [50] to reproduce this orphan *γ*-ray flare by introducing an extreme hardening of the electron spectrum, along with a location of the emission region much closer to the BH and accretion disk. A hard electron spectrum *ne*(*γ*) ∝ *γ*−1 up to a cut-off energy of a few 1000 is invoked, which may be difficult, but not impossible, to achieve with standard particle acceleration mechanisms. Asano and Hayashida [51] employ a time-dependent one-zone model with second-order Fermi acceleration, where an enhanced acceleration efficiency leads to a hardening of the electron spectrum, and a significant reduction of the magnetic field is required to suppress a simultaneous optical flare. While their model represents the *γ*-ray spectrum during the flare well, it does predict a non-negligible optical synchrotron flare accompanying the *γ*-ray flare. A similar strategy, based on an analytical solution to the steady-state electron distribution, was adopted by Lewis et al. [52], also requiring a significant reduction of the magnetic field to suppress a simultaneous optical synchrotron flare. Yan et al. [53] modelled the orphan-flare SED using a time-dependent single-zone model with rapid electron cooling. However, it is unclear whether a transition from the quiescent to this flaring state may

be produced in a natural way. Lepto-hadronic models naturally de-couple the (protoninitiated) high-energy emission from the (electron-initiated) synchrotron radiation and therefore offer an alternative way of reproducing orphan *γ*-ray flares. Paliya et al. [54] used the time-dependent lepto-hadronic model of Diltz et al. [55] to model the December 2013 orphan *γ*-ray flare of 3C279. They also considered the possibility of a two-zone model, with a small emission region emitting an SED with very large Compton dominance, emerging during the time of the orphan *γ*-ray flare.

A representative simulation of the shock-in-jet synchrotron mirror scenario was then used for a generic study of the expected spectral variability patterns. X-ray and *γ*-ray light curves as well as radio and optical light curves are significantly correlated with each other, while radio and optical light curves are significantly anti-correlated, with radio and optical (synchrotron) dips accompanying the high-energy flare resulting from more efficient radiative cooling of the electrons. The response in the radio light curves is found to be delayed by ∼15 h with respect to other bands. In the scenario investigated here, where no changes to the diffusive shock acceleration process along the entire evolution of the flare are assumed, significant spectral hysteresis is not expected, but a mild harder-when-brighter trend in most wavebands is found.

While it is found that the specific December 2013 orphan *γ*-ray flare of 3C279 can not be successfully reproduced with this scenario, it may be applicable to other, more moderate orphan *γ*-ray flares. Especially the expected anti-correlation between Compton and synchrotron wavebands may serve as a smoking-gun signature of this scenario. The failure of the shock-in-jet synchrotron mirror model for the December 2013 flare of 3C279 was primarily caused by the fact that the shock-accelerated, relativistic electrons were already in the fast-cooling regime and radiating very efficiently, as was required by the fit to the quiescent state of 3C279. If the quiescent emission of a blazar can be produced by a less radiatively efficient configuration, then the increase in radiative efficiency in the synchrotron mirror scenario may lead to substantially higher-amplitude flares. A systematic study of different scenarios and applications to other sources will be presented in a forthcoming publication.

**Funding:** The work of M.B. is supported by the South African Researc Chairs Initiative of the National Research Foundation (any opinion, finding, and conclusion or recommendation expressed in this material is that of the authors, and the NRF does not accept any liability in this regard) and the Department of Science and Innovation of South Africa through SARChI grant no. 64789.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


*3.2. DNA Oxidation*

DNA oxidation was assessed by quantitative immunostaining of the oxidized guanine in the nucleus (8-OHdG). Immunofluorescent staining was detected in both epithelial and interstitial cells of the different regions of the epididymis (initial segment, caput, corpus, and cauda) (Figure 2). In 18-month-old mice, staining was more intense in the cauda epididymidis than in the rest of the tissue; further, this increase was more pronounced in the distal cauda epididymidis than in the proximal cauda epididymidis. Quantification of nuclear 8-OHdG immunofluorescent staining (Figure 2 bar graph) indicated that the major increase in DNA oxidation was observed in the cauda epididymidis of 18-month-old *Sod1*−/− mice. Statistical analysis by 3-way ANOVA (age, genotype, region, Table 2) revealed that the main source of variation in the DNA oxidation levels was the *Sod1*−/− genotype; in addition, there were interactions between epididymal regions and age and between age and genotype.

\* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001, \*\*\*\* *p* < 0.0001.

**Figure 2.** Aging induces an accumulation the nucleic acid damage in the epididymis (initial segment, caput, corpus, and cauda epididymides) which is enhanced in *Sod1*−/−mice. Representative pictures of sections of the epididymis from 3-month-old and 18-month-old wild-type and *Sod1*−/−mice after immunostaining of oxidized nucleic acids (8-OHG, red); the nuclei are counterstained with DAPI (blue). Immunostaining negative controls (no primary antibody) are displayed for the caput and cauda epididymides. Scale bar: 40 μm. 8-OHG staining was quantified separately in the epididymal epithelium (clear histograms) and in the interstitial tissue (dashed histograms). \*\* *p* < 0.01, \*\*\* *p* < 0.001 (3-way ANOVA, n = 4–5).