2.4 Steroid receptors and their ligands
2.4.1 Androgen receptor in prostate cancer
The androgen receptor (AR) is a member of the NR superfamily, expressed in a diverse range on tissues (Davey and Grossmann, 2016). Its natural ligand
dihydrotestosterone (DHT) is a steroid hormone that regulates the development of male reproductive system and the secondary sexual characteristics.
The AR has a similar overall structure and domain structure as the other NRs.
However, AR has an exceptional C-terminal extension that is found to be folded on the surface of the LBD when the helix 12 adopts the agonistic conformation. Like many steroid hormone receptors, AR functions as a homodimer. The C-terminal extension blocks the formation of the type of dimer that is seen in the crystal structures of RXR heterodimers by folding onto the dimerization surface (Bledshoe et al. 2002). The structure of the homodimer is different from many other
structures, e.g PPAR heterodimers with RXR, but it does resemble the structure of the GR homodimer (Nadal et al., 2017).
Prostate cancer is common in men over 50 years of age. Genetic alterations commonly occur in signalling pathways that promote growth as well as those blocking apoptosis in prostate cancer cells. Some of these pathways are regulated by AR, and some are able to alter the structure of the AR itself. This, along with the fact that the activity of AR can be regulated with small molecules, makes the AR an important target for the treatment of prostate cancer. (Shafi et al., 2013)
Prostate cancer is often treated with androgen deprivation therapy (ADT), which means deprivation of the AR of testosterone. This is achieved by surgical or medical orchidectomy. AR antagonists can be used in addition to achieve complete androgen blockade. (Davey and Grossmann, 2016)
After an initial response, most tumors become refractory to ADT. This state is called hormone refractory or castration resistant prostate cancer (CRPC). There are multiple mechanisms leading to the loss of response, including an increase in AR expression, AR point mutations and other variants, AR function modulation by changes in cell signalling pathways, and changes in coregulator proteins or androgen biosynthesis. (Fujita and Nonomura, 2019; Shafi et al., 2013) 2.4.2 Antiandrogens and AR mutations
Compounds that compete with androgens for binding to AR and inhibit its activation are called antiadrogens. Steroidal antiandrogens like cyproterone acetate were used to treat advanced prostate cancer because of their ability to bind to AR and block DHT and testosterone binding (Figure 8). However, steroidal antiandrogens had side effects which limited their usability. Steroidal
antiandrogens were commonly replaced by nonsteroidal compounds in prostate cancer therapy because of their more favourable safety profile. First generation nonsteroidal antiandrogens flutamide and nilutamide and second generation
compound bicalutamide all have similar, relatively low affinities for the AR (Figure 9). This means that DHT can still activate AR to stimulate the growth of prostate cancer cells. (Crawford et al., 2018)
Figure 8. Structures of steroidal AR agonists dihydrotestosterone and R1881 and antagonist cyproterone acetate.
Enzalutamide (Figure 9) is a third generation antiandrogen which, compared with bicalutamide, has five to eight fold higher binding affinity for the AR. Apalutamide has a similar structure as enzalutamide. Darolutamide is structurally different than the other two third generation antiandrogens. When compared to other two, its advantage is its lower penetration across the blood–brain barrier. It is also active against some of the common AR mutants, including the F876L mutation that makes the tumor resistant to enzalutamide and apalutamide. (Crawford et al., 2018; Tran et al., 2009)
Figure 9. Structures nosteroidal of AR antagonists.
The use of antiandrogen commonly leads to resistance, by mechanisms including AR overexpression and AR mutations. Point mutations are found in approximately 15–30 % of patients with CRPC, most commonly in the ligand binding pocket (Fujita and Nonomura, 2019). For example, many patients respond well at first to
treatment with enzalutamide, but later they develop resistance. One mechanism for the resistance is the appearance of the F876L mutation in the AR ligand binding pocket, which allows the AR to use enzalutamide as an agonist to drive tumour growth. However, the AR with the F876L mutation still responds to bicalutamide. In contrast, the W741C mutation turns bicalutamide into an agonist, but this is not the case for enzalutamide. Mutation T877A confers AR resistance to
hydroxyflutamide (Figure 9) and cyproterone acetate (Figure 8), and AR with T878A can be activated with flutamide, bicalutamide and exalutamide, but also with progesterone and estrogen. Because of these kinds of differences in activity profiles of antiandrogens in AR mutants, it has been suggested that the use of structurally distinct compounds in combination or in series might represent an
efficient strategy for combating the resistance caused by AR mutations. (Fujita and Nonomura, 2019; Romanel et al., 2015; Bambury and Scher, 2015)
2.4.3 Rational design of AR ligands
Because the resistance to a certain structure is a common problem in the use of antiandrogens, much work has been done in developing new types of
antiandrogens. The resistance to the antiandrogens flutamide and bicalutamide was associated with over-expression of the AR, suggesting that a more potent AR inhibitor without partial agonist activity could be used in these cases. This was the starting point for the development of enzalutamide. The new structure was found by developing a library of compounds based on a known AR binder and testing them in prostate cancer cell lines. (Bhattacharya et al., 2015)
Another approach was taken by Bassetto et al. Their work started with the observation that one of the most common mutants of the AR, W741L, converted bicalutamide into an agonist. The reason for this phenomenon lies in the flexible structure of bicalutamide. Mutation of W741 to a smaller residue creates more space within the ligand binding pocket of AR, and bicalutamide can be twisted to fill that space, thus leading to the correct positioning of the AR helix 12 for activation. Bassetto and coworkers aimed to design a structure that would resemble bicalutamide, but would be more rigid to avoid the loss of antagonistic activity caused by mutations. For this purpose, the structural features of
bicalutamide were combined with another, more rigid, AR binder to create a new compound (Figure 10). (Bassetto et al., 2017)
Figure 10. Rational design of new antiandrogens. Figure redrawn from (Bassetto
Prekovic et al. used docking to examine the effects of AR mutations on ligand responses (Prekovic et al., 2016). For example, the F877L mutation expands the ligand binding pocket of AR and allows enzalutamide to bind and behave like an agonist. Similar conclusions had been made earlier based on MD simulations (Balbas et al., 2013). The authors suggested that these observations should be taken into account when designing the next generation antiandrogens, which could maintain their antagonist activity even in mutated receptors.
2.4.4 The estrogen receptor
The estrogen receptor (ER) is one of the steroid receptors belonging to the NR superfamily. It has two subtypes, named ERα and ERβ. The subtypes are highly homologous in the DBD, but only 56 % in the LBD. In the ligand binding pocket, however, there are only two amino acids that are different between ERα and ERβ.
Human ERα has 595 amino acids (66 kDa) and ERβ contains 530 amino acids (54 kDa). ERs have a C-terminal F-domain (see Figure 3), which affects the activity of the ERs and may contribute to the selective activation of specific target genes. The natural ligand for both isotypes is 17β-estradiol (E2). ERs function as homodimers.
(Liu, J., 2020; Begam et al., 2017; Moon et al., 2009; Dahlman-Wright et al., 2006) In some organs, both isotypes are expressed at similar levels, but ERα or ERβ may predominate in some tissues or different cell types in the same tissue. ERα is mainly expressed in uterus, testes, bone, breast, liver and white adipose tissue while ERβ tends to be expressed in colon, testes and vascular endothelium. In ovary, prostate and brain, both isotypes are expressed but in different regions. ERs have a similar domain structure as the other receptors in the NR superfamily.
(Begam et al., 2017; Moon et al., 2009)
ERα and ERβ have their own roles and functions in the immune, skeletal, cardiovascular and central nervous systems. Their activation can even have opposite effects, and the final outcome depends on the balance between ERα and ERβ signalling. Despite the high similarity of the ligand binding pockets, it has been possible to develop ligands that are selective for either ERα or ERβ. (Begam et al., 2017; Moon et al., 2009)
2.4.5 Selective estrogen receptor modulators
Compounds which have high affinity for ERs but not for other steroid hormone receptors, are called selective estrogen receptor modulators, or SERMs. The interesting feature of these compounds is that they can display tissue specific agonistic or antagonistic activity, that is, behave as agonists or antagonists
depending on the tissue type. This kind of feature is of major interest in drug development because it offers a possibility to avoid undesired side effects outside of the target tissue. A good example is tamoxifen, which acts as an antagonist in the breast but as an agonist in bone, so it can be used simultaneously to treat osteoporosis and prevent breast cancer. These opposite outcomes are a result from the cell surroundings in different tissues, varying levels of expression of coactivator and corepressor proteins and phosphorylation of ER. A ligand binding to ER can change the shape of ER and in this way, also the interactions it can have with its coregulators. However, e.g. high level of expression of coactivators may convert the equilibrium from antiestrogenic behaviour in one cell type to
estrogenic behaviour in another type of cells. (Maximov et al., 2013; Dutertre and Smith, 2000)
Docking and MD have been used in attempts to predict the agonistic or
antagonistic properties ER binders. (Cotterill et al., 2019; Puranik et al., 2019; Ng et al., 2014) Docking into the agonistic and antagonistic crystal structures of ER can be used to distinguish agonistic and antagonistic binders. However, it is unlikely that these methods can accurately predict the very subtle differences in agonistic and antagonistic behavior, especially in the case of tissue specific effects.
Figure 11. Structures of estradiol and selective estrogen receptor modulators. All
2.4.6 Clinical use of ER antagonists
The first generation SERMs, tamoxifen and toremifene (Figure 11), are still in clinical use in the treatment of breast cancer. In addition to their antagonist activity in breast and agonistic activity in bone, tamoxifene also has agonistic activity in uterus; with long term use, it may increase the risk of endometrium cancer. Toremifene resembles tamoxifene in terms of both structure and
biological activity. A second generation compound, raloxifene, is indicated for the treatment and prevention of osteoporosis. It does not have a similar agonistic effect in uterus i.e. it differs from tamoxifen and toremifene. (Liu, J., 2020; Dutertre and Smith, 2000)
Basedoxifene is a third generation SERM which has a strong antagonistic effecs in breast and endometrium. It is used to treat vasomotor symptoms in
menopausal women. Ospemifene is structurally similar to tamoxifen and toremifene with similar biological properties. Lasofoxifene, another third generation SERM, is used in the treatment of osteoporosis and it is being
evaluated to treat breast cancer. Arzoxifene failed in phase III trials and has been discontinued. (Li, H. et al., 2020; Liu, J., 2020; Dutertre and Smith, 2000)
3 CO PUTATIONAL ETHODS
3.1 MD SIMULATIONS
3.1.1 Basic concepts
Structural information of biomolecules from crystallographic or nuclear magnetic resonance (NMR) measurements is needed for MD simulations. When the
positions of the atoms in the system are known, it is possible to calculate the forces exerted on each atom by other atoms in the system. The forces are used to calculate the new positions of the atoms after a short period of time, and the new forces exerted on the atoms. By repeating this sequence of calculations, a
trajectory of atoms moving over time is produced. Thus, it is possible to calculate the movements of the molecule. The bonds between atoms cannot form or break, so it is not possible to model chemical reactions with a classical MD simulation.
(Hollingsworth and Dror, 2018; Durrant and McCammon, 2011)
The forces described above are calculated using a molecular mechanical force field. The term ‘force field’ refers to a mathematical equation and its parameters, which are used to calculate the potential energy of the system. Typically, the force field contains definitions for different atom types, parameters for bonds, angles and torsions as well as for non-bonded interactions. An example of a force field equation for calculating the forces in the system is shown in equation (1).
(Monticelli and Tieleman, 2013)
𝐸𝑡𝑜𝑡𝑎𝑙= ∑ 𝑘𝑑
The first term in the equation represents the bond stretching with a harmonic potential. kd represents the force constant and d0 the equilibrium bond length.
Angle bending and improper dihedrals are described analogously. The improper dihedrals term is used for out-of-plane bending of planar groups.
Dihedral-term is applied on torsion angles for atoms separated by three bonds.
In this term, k∅ is a force constant related to the barrier height. The number of
minima in the energy function is determined by n. The phase factor ∅0 determines the positions of the minima.
The last two terms represent non-bonded interactions. The non-bonded interactions are usually calculated for atoms which are not connected with bonds or which are separated by at least three bonds, but a scaling factor can be applied on non-bonded interactions, which are also covered by torsion terms. The second-last term is called Lennard-Jones 12-6 potential, and it represents the van der Waals component of the potential energy. In this term, σ is the distance defined for each pair of atom types, in which the potential between the two atoms in question is zero. In shorter distances there is a repulsive force, and in longer distances there is an attractive force, the strength of which is determined by ε. The last term in the equation represents the electrostatic interactions between atoms with (partial) charges (q). ϵ is the dielectric constant.
The force field parameters are developed to reproduce quantum mechanical or experimentally observed geometries and energies in chosen simulation
conditions. A number of force fields are currently available for simulation of biological macromolecules, e.g. Amber (Cornell et al., 1995), Gromos (Oostenbrink et al., 2005), OPLS (Harder et al., 2016) and CHARMM (Mackerell, 2004). The major force fields applied in MD simulations have also parameters for various polar and apolar solvents, common organic molecules, proteins, carbohydrates, nucleic acids and some common phospholipids. (Monticelli and Tieleman, 2013)
The force field is used to calculate forces on each atom, and the Newtonian equation of motion is then used to calculate the accelerations and new positions of the atoms after a short time step. Verlet (Verlet, 1967) and leap-frog (Hockney, 1970) algorithms are common methods to numerically integrate the equation in a MD simulation. (Allen, 2004; Leach, 2001)
The MD simulations are usually performed in a certain temperature and pressure to mimic in vivo or in vitro conditions. The temperature in the MD simulation is calculated from the kinetic energy of the system. Controlling the temperature of the system can be done by scaling the velocities of the atoms in the system. Commonly used algorithms for temperature control in MD simulations are Berendsen (Berendsen et al., 1984) and Nosé-Hoover (Hoover, 1985; Nosé, 1984) algorithms. Constant pressure conditions are maintained by changing the volume of the system. Berendsen (Berendsen et al., 1984) and Parrinello–Rahman (Parrinello and Rahman, 1981; Parrinello and Rahman, 1980) are commonly used algorithms for pressure control. (Leach, 2001)
3.1.2 Challenges and limitations of MD simulations
MD simulations are used to study dynamic properties of molecules. This kind of knowledge is difficult to get with other methods. However, a series of calculations for a system with many atoms requires a lot of computational capacity, and this task is currently not possible without a simplification of the description of the system. MD simulation uses a force field to approximate the quantum mechanical wavefunction, reducing the computational task from quantum mechanics to classical mechanics, ignoring the electronic motions. (Vanommeslaeghe et al., 2014; Guvench and MacKerell, 2008; Leach, 2001)
The outcome of the MD simulation is crucially dependent on the success of the parametrization of the force fields, which are constantly developed futher for better accuracy. Different force fields may have a different approach on their parametrization and the equations may contain different terms. The performances of the force fields may also differ, when the geometries they produce are
compared to e.g. NMR measurements. A force field parametrized for a certain specific purpose will perform better in that purpose, compared to a more general-purpose force field. In addition, a classical MD force field uses partial charges assigned on each atom before the simulation starts (Jorgensen, 2007).
Ignoring electronic motion means that electronic polarization on response to surroundings of the atoms is not captured. Polarizable force fields have been developed which improve accuracy, but with the cost of increasing computational time (Cieplak et al., 2009). (Rahman et al., 2020; Henriques et al., 2015; Beauchamp et al., 2012; Durrant and McCammon, 2011)
Interesting biological processes, such as ligand binding and conformational changes of protein, happen on timescales from nanosecond to microsecond or longer. On the other hand, the length of the time steps in a classical all-atom MD simulation must be short, typically 1–2 femtoseconds, to avoid instability, so a huge number of time steps are required for a long timescale MD simulation. Thus, the computational capacity currently available also limits the length of the
simulations. (Hollingsworth and Dror, 2018; Kmiecik et al., 2018)
3.2 FREE ENERGY CALCULATIONS
Free energy calculations can be used to provide estimates of free energy
differences between two states of a system, for example, a protein with a bound and unbound ligand. Binding of a drug to its target protein is required for its therapeutic effects, and thus it is important in drug design to be able predict
binding affinities accurately and reliably. There are several ways to calculate the free energy on binding. Molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) (Poli et al., 2020; Wang et al., 2017; Srinivasan et al., 1998) and thermodynamic integration (TI) (Kollman, Peter, 1993) methods are shortly described here. (Shirts et al., 2010)
MM-PBSA is a widely used method to calculate the free energy of binding (ΔGbind) of a small ligand to a protein in aqueous solvent. It is an end-point method, which requires simulations of only the bound and unbound states. The Gibbs free energy of binding is given by:
ΔG
bind= ΔH -TΔS
(2)where ΔH is the difference in enthalpy between the bound and unbound states and -TΔS is the entropy difference in a temperature T. In MM-PBSA method, the equation (2) takes the form:
ΔG
bind= ΔE
MM+ ΔG
solv-TΔS
(3)Th
e gas-phase interaction energy EMM between the ligand and the protein is calculated from the molecular mechanical force field. Gsolv is the solvation free energy which can be calculated as a sum of its polar and non-polar components.The polar contribution is calculated using Poisson–Bolzmann implicit solvent model (Honig and Nicholls, 1995) and the non-polar component based on solvent accessible surface area (Srinivasan et al., 1998). Generalized Born -method can be also used for the polar component, in which case the method is called MM-GBSA (Genheden and Ryde, 2015). The term -TΔS in equation (3) is the entropy change of the solute associated with the ligand binding. The entropic term is
computationally expensive to calculate, and it is often not included if the aim is not to calculate absolute binding free energies but rather, to compare the free
energies of binding of series of ligands. In practice, a free energy calculation with MM-PBSA approach utilizes MD simulations, which are used to get the snapshot structures for calculation of the energy terms, which are calculated as an average over the snapshots. (Poli et al., 2020; Wang et al., 2017; Shirts et al., 2010)
TI method (Kollman, Peter, 1993) is used to calculate the relative free energy difference of two states, for example, the difference of binding free energies of two ligands. The system is changed from one state (λ=0) to another state (λ=1).
The coupling parameter λ is used to link the two physical states via a pathway of