# Location Aided Energy Balancing Strategy in Green Cellular Networks

Jianjun Yang  Department of Computer Science
University of North Georgia
Email<EMAIL_ADDRESS>Bryson Payne  Department of Computer Science
University of North Georgia, USA
Email<EMAIL_ADDRESS>Markus Hitz  Department of Computer Science
University of North Georgia, USA
Email<EMAIL_ADDRESS>Zongming Fei  Lab for Advanced Networking
Department of Computer Science
University of Kentucky
Email<EMAIL_ADDRESS>Le Li  David R.heriton School
of Computer Science
University of Waterloo, Canada
Email<EMAIL_ADDRESS>Tongquan Wei  School of Information Science
and Technology
East China Normal University, China
Email<EMAIL_ADDRESS>
###### Abstract

Most cellular network communication strategies are focused on data traffic
scenarios rather than energy balance and efficient utilization. Thus mobile
users(cell phones) in hot cells may suffer from low throughput due to energy
loading imbalance problem. In state-of-art cellular network technologies,
relay stations extend cell coverage and enhance signal strength for mobile
users. However, busy traffic makes the relay stations in hot area run out of
energy quickly. In this paper, we propose an energy balancing strategy in
which the mobile nodes are able to dynamically select and hand over to the
relay station with the highest potential energy capacity to resume
communication. Key to the strategy is that each relay station merely maintains
two parameters that contains the trend of its previous energy consumption and
then predicts its future quantity of energy, which is defined as the relay
station’s potential energy capacity. Then each mobile node can select the
relay station with the highest potential energy capacity. Simulations
demonstrate that our approach significantly increase the aggregate throughput
and the average life time of relay stations in cellular network environment.

##  I Introduction

During the past several years, there has been tremendous growth in cellular
networks. With the introduction of Android and iPhone devices, use of ebook
readers such as iPad and Kindle and the success of social networking giants
such as Facebook, WeChat and QQ, the number of customers and the demand for
cellular traffic have escalated significantly[26][22] [27][30]. The huge
amount of customers and very high volume data transmission result in serious
problems of energy consumption[2].

The rising energy costs of operating cellular networks have led to an emerging
trend of addressing energy efficiency and energy balance utilization amongst
the network operators and regulatory bodies such as 3GPP and ITU[3]. This
trend has stimulated the interest of researchers in an innovative new research
area called “green cellular networks”. In this regard, the European Commission
has recently started new projects to address the energy issue of mobile
communication systems, such as “Energy Aware Radio and Network Technologies
(EARTH)”, “Towards Real Energy-efficient Network Design” and “Cognitive Radio
and Cooperative strategies for Power Saving in Multi-standard Wireless
Devices” [4]. Energy in cellular network is a vast research discipline that
needs to cover all the layers of the protocol stack and various system
architectures and it is important to identify the fundamental trade-offs
linked with energy efficiency and the overall performance [5].

Scholars have addressed four key issues in terms of energy efficiency with
network performance [3][23] . They are deployment efficiency (balancing
deployment cost, throughput), spectrum efficiency (balancing achievable rate),
bandwidth (balancing the bandwidth utilized) and delay (balancing average end-
to-end service delay). To address the challenge of increasing power efficiency
in future cellular networks and thereby to maintain profitability, it is
crucial to consider various paradigm-shifting technologies, such as energy
efficient wireless architectures and protocols, efficient base station(BS)
redesign, smart grids, opportunistic network access or cognitive radio,
cooperative relaying and heterogeneous network deployment based on smaller
cells.

In addition to academia, governments and industries have recently shown keen
concerns on the critical issues related to energy efficiency and its balance
utilization and in the ICT(Information and Communication Technology) area.
However, as studied in the recent literature [6][7], most of the techniques
applied to current mobile networks have been designed by taking into account
non-energy-related factors, such as throughput, Quality of Service (QoS),
availability, scalability, and so on[34][29][35][32][21][33][28]. Moreover, in
real-world systems, mobile users (MUs) are not evenly distributed across
cells, resulting in that MUs in a hot cell will be affected by the load
imbalance and they might unable to get services. Meanwhile, as novel network
architectures that include picocells, hierarchical cells and femtocells
emerge, the density of base stations and mobile users is becoming larger, and
the cells go smaller. The appearance of this high density cellular networks
introduces more variety of load across different cells and makes the load
imbalance problem more serious[36][37].

The new generation cellular networks allow mobile users to connect relay
station and then connect to base station. The busy traffic near a relay
station often makes the energy of the relay station go down very quickly.
However, the relay station with low traffic is idle with high energy. In order
to balance the energy utilization among different cells, it is needed to
transfer the over-loaded traffic from hot cells to neighboring cooler ones.
The challenge is how to balance the energy utilization of the relay stations
in order to get the best trade-off among all the relay stations in a cellular
network.

In this paper, we develop a statistical parameter based energy balance
utilization algorithm, in which each relay station maintains the acceleration
and the variance of energy consumption that represent the station’s historical
energy consumption acceleration and the variance. Together with its current
energy quantity, the relay station is able to predict its future energy
quantity, which is considered as its potential energy capacity. Before payload
data transmission between a mobile user and a relay station, the mobile user
disseminates a message to a portion of its neighbor relay stations depending
on their geographic locations . Then it selects and hands over to the relay
station with highest potential energy. Because each relay station only needs
to maintain two parameters, the overhead of this scheme is very low. Our
simulation results illustrate that our approach significantly increases the
aggregate throughput in the network and the average life time of the relay
stations compared with existing approaches.

The rest of the paper is organized as follows. Section II discusses the
related research on this topic. Section III proposes a novel method that
select the best relay station. We evaluate the proposed schemes by simulations
and describe the performance results in Section IV. Section V concludes the
paper.

##  II Related Work

Various energy utilization strategies in cellular networks were proposed and
scholars developed various ways to solve energy utilization
questions[38][25][39][40]. The schemes are classified into several categories.
The first one contains the strategies based on channel borrowing from cooler
cells[10]. The second category includes the strategies based on BS
selection[11]. The third category contains strategies based on power control
and cell breathing[11][12]. The last type consists of strategies on relay-
assisted traffic transfer[13]. The basic idea of channel borrowing is to
borrow a set of channels from cooler cells to hot cells. However, this will
change the pre-defined spectrum reuse pattern and introduce more co-channel
interference[8]. There are also a great deal of research combining multi-hop
wireless network (MWN) with infrastructure wireless networks, forming
architecture of multi-hop cellular network (MCN). In MCNs, relay stations
(RSs) are network components that are dedicated to storing and forwarding data
received from BSs to MUs, and vice versa. Deploying relay stations can extend
the coverage and enhance the signal strength, which clearly help improve the
performance for MUs near the edge of the cell. The appearance of MCNs also
provides another method to solve the load imbalance among different cells.
That is to transfer over-loaded traffic from hot cells to cooler cells by RSs.
Compared with previous discussed dynamic load balancing schemes, relay-based
load balancing schemes are more flexible and will introduce less interference.
In[6], a mobile-assisted call admission scheme is proposed to achieve load
balancing in cellular networks, which require an ad-hoc overlay network on the
cellular network. The authors divided the channels into two parts, one for the
ad-hoc overlay network, the other for the cellular network. In[7][8], the
authors proposed dynamic load balancing schemes in the integrated cellular and
ad-hoc relaying systems (iCAR). The ad-hoc relaying stations (ARS) compose an
overlay ad-hoc network, which can help relay traffic among different cells. As
in iCAR systems, ARSs work on ISM-band channels, it is pointed out that the
performance of iCAR systems will depend on the number of available ISM-band
channels, and interference in ISM-band could affect the performance of dynamic
load balancing.

A. Alam et al. [15] investigates dynamic traffic-aware BS switching modes,
where the BS can alter its operating modes between standard BS operations and
switching to relay station (RS) mode, the so-called BS-RS Switching model.
Depending on the traffic fluctuations, load profiles are divided into two
categories. One is Zero-to-medium traffic period when a BS switches to the RS
mode and turns off all its high-power consuming equipment. The second is Peak
traffic period when all BS are fully active. The rationale for switching from
a BS to RS mode is to ensure those MS that would be served by the switched off
cell and may suffer deep fading, are still be able to receive the same QoS.
Furthermore, since the propagation distance has been shortened between the MS
and serving BS via the back haul connection, the required MS transmit power is
concomitantly reduced compared with the BS sleep and cell zooming technique
[5].

K. Xu and M. Zhou [16][17] developed a realistic and representative energy
model based on RF transceiver. As a typical IEEE 802.15.4 RF transceiver mode,
each wireless node is equipped with multiple discrete transmit power levels.
Based on this energy model, they developed a power control policy, which
configure the transmit power level as a function of transmission distance such
that the energy cost is minimized. Then, the optimized energy balanced chain
model is proposed to determine the optimal traffic flow distribution.

##  III Problem Formulation

###  III-A The Basic Idea

In cellular networks, each base station covers a number of cell phones in a
hexagonal cell. We consider a cellular network containing a set of base
stations(BS), relay stations(RS) and mobile users(cell phones)(MU), denoted by
$\left\\{B_{1},B_{2},...,B_{s}\right\\}$,
$\left\\{R_{1},R_{2},...,B_{t}\right\\}$, and
$\left\\{m_{1},m_{2},...,m_{p}\right\\}$. Each BS is located at the center of
the cell. Each RS is located at a boundary of cells. The traditional
communication between two mobile users $m_{i}$ to $m_{j}$ is that $m_{i}$
communicates to its base station $B_{i}$, then $B_{i}$ communicates to $B_{j}$
which covers $m_{j}$, and then $B_{j}$ communicates to $m_{j}$. In new
generation cellular network, the relay stations are used to assist
communications between cell phone to cell phone or cell phone to base station.
Fig. 1 shows an instance with several base stations, relay stations and cell
phones, where the central cell is a hot cell with a large number of cell
phones and busy data traffic and the other cells are cooler cells with less
cell phones. Note that the cell phones in all figures of this paper represent
active cell phones.

Figure 1: A topology of a cellular network.

Data transmission takes place when a cell phone(initiator) intends to talk to
another cellphone(receiver). The scenario is simple if the two are in the same
cell[9]. We consider the circumstance that the receiver is in another cell. In
Fig. 1, the communication between $m_{1}$ and $m_{2}$ is such an instance, in
which data is relayed from $m_{1}$ to $R_{1}$, then from $R_{1}$ to $B_{1}$,
and continues until it reaches $m_{2}$. In the scenario, $R_{1}$ is hot with a
lot of cell phones nearby. Hence its energy is going down very quickly.
However, $R_{2}$ covers less cell phones and its energy is going down slowly.
Here comes up with a question, can $m_{1}$ find $R_{2}$ or even $R_{3}$ to
relay the message to bypass the hot relay station $R_{1}$?

Since a mobile cell phone is able to select a relay station and hand over to
it[9], we develop an approach to let each mobile user find the RS with highest
potential energy and hand over to it. We start our basic idea by a simplified
case as Fig. 2. Suppose $m_{1}$ intends to talk to $m_{2}$. If $m_{1}$ selects
$R_{1}$ to conduct relay, it may not be a good choice because $R_{1}$ is busy
and may run out of energy shortly. $R_{1}$, $R_{2}$ and $R_{3}$ are all
possible relay stations for $m_{1}$. Suppose the residual energy of $R_{1}$ is
900 $kWh$, $R_{2}$ is 800 $kWh$ and $R_{3}$ is 850 $kWh$ at current moment. We
further assume that the historical residual energy vector up to this moment of
$R_{1}$ is […, 2000, 1500, 1300, 900], $R_{2}$ is […, 900, 870, 830, 800], and
$R_{3}$ is […, 1400, 1200, 1000, 850]. A simple strategy will let $m_{1}$
select the one that has the highest energy at the moment. Seemingly, $R_{1}$
is the best choice. However, since $R_{1}$ is in a hot area serving a large
number of mobile users and its energy is going down dramatically, while
$R_{2}$ are serving much less mobile users and its energy goes down slowly,
$m_{1}$ should select $R_{2}$.

Our goal is to let $m_{1}$ select the relay station that will have the highest
potential energy. In our approach, each RS can predict its future energy right
after this time slot based on its historical and current energy, which is
defined as potential energy of the RS in this paper. In this example, our
strategy will let $m_{1}$ choose $R_{2}$.

Figure 2: The scenario to illustrate the main idea.

###  III-B Computation of Potential Energy

The mobile initiator needs to know which relay station is the one that has the
potential highest energy. The energy varying trend is essential to determine
an RS’s potential energy. Suppose the current time is $k$. A naive method on
predicting the energy at time $k+1$ is as follows. Let RS $R_{i}$ save all its
historical and current (time $k$) energy $E_{i_{0}}$, $E_{i_{1}}$,
…,$E_{i_{k}}$ in a vector. Then it can predict its energy for time $k+1$ by
curve-fitting in numerical analysis that approximates its moving trend. The
limitations of this method are to spend too much space to save the energy
values in the vector and too complex for curve-fitting computation.

We propose a statistical based strategy to let each relay station predict its
future energy. In this scheme, each RS only needs to maintain two parameters.
One is its potential acceleration and the second is the variance of the
acceleration. The two parameters store how its residual energy changes.
Together with the RS’s current energy, its future energy capacity can be
calculated and is considered its potential energy. Thus, mobile user $m_{1}$
is able to figure out the best RS subsequently.

We define the following notations to represent the terms regarding our
approach for node $R_{i}$.

$E_{i_{k}}$: The measured energy at time $k$.
$\hat{E_{i_{k}}}$: The energy predicted at time $k$.
$a_{i_{k}}$: The acceleration measured at time $k$. It indicates how energy
changes during a time slot.
$\hat{a_{i_{k}}^{-}}$: The acceleration at time $k$ evolved from time $k-1$.
$\hat{a_{i_{k}}}$: The potential acceleration at time $k$.
$v_{i_{k}}$: The variance of acceleration updated at time $k$.
$v_{i_{k}}^{-}$: The variance of acceleration at time $k$ evolved from time
$k-1$.
$\epsilon$: The error or noise in the process.
$B_{i_{k}}$: The blending factor at time $k$.

At time $k$, $R_{i}$ measures its energy $E_{i_{k}}$. And then it computes
$a_{i_{k}}$ as its measured acceleration by

$a_{i_{k}}=(E_{i_{k}}-E_{i_{k-1}})/\Delta t~{}~{}~{}(1)$

$R_{i}$ updates $\hat{a_{i_{k}}^{-}}$ and $v_{i_{k}}^{-}$ in order to keep its
historical energy to predict its future energy.

$\hat{a_{i_{k}}^{-}}=\hat{a_{i_{k-1}}}~{}~{}~{}(2)$
$v_{i_{k}}^{-}=v_{{i_{k}-1}}+\epsilon~{}~{}~{}(3)$

$R_{i}$ also computes the blending factor $B_{i_{k}}$, which indicates how
much the acceleration changes from last time to current time.

$B_{i_{k}}=v_{i_{k}}^{-}(v_{i_{k}}^{-}+\epsilon)^{-1}=v_{i_{k}}^{-}/(v_{i_{k}}^{-}+\epsilon)~{}~{}~{}(4)$

Once $R_{i}$ obtains the blending factor $B_{i_{k}}$ and the evolved
acceleration $\hat{a_{i_{k}}^{-}}$, it knows how much the acceleration changes
and the evolved acceleration. Additionally, $R_{i}$ considers the measured
acceleration $a_{i_{k}}$. Then it calculates its potential acceleration
$\hat{a_{i_{k}}}$. This acceleration will be used to predict its energy of
time $k+1$.

$\hat{a_{i_{k}}}=\hat{a_{i_{k}}^{-}}+B_{i_{k}}(a_{i_{k}}-\hat{a_{i_{k}}^{-}})~{}~{}~{}(5)$

$R_{i}$ updates the variance of acceleration for future utilization.

$v_{i_{k}}=(1-B_{i_{k}})v_{i_{k}}^{-}~{}~{}~{}(6)$

Finally, the energy predicted of time $k+1$ is

$E_{i_{k+1}}=E_{i_{k}}+\hat{a_{i_{k}}}\Delta t~{}~{}~{}(7)$

At a certain time $k$, $R_{i}$ only needs to measure its energy $E_{i_{k}}$
and record two parameters $\hat{a_{i_{k-1}}^{-}}$ and $v_{i_{k-1}}$. Then it
can predict its energy at time $k+1$ by the calculation, which is the
potential energy capacity of $R_{i}$.

###  III-C Algorithm Description

In our approach, when a user $m_{1}$ intends to communicate to $m_{2}$, it
selects the relay station with highest potential energy. Apparently, $m_{1}$
does not need to consider the relay stations in the opposite direction $m_{1}$
$\rightarrow$ $m_{2}$, such as $R_{4}$, $R_{5}$ and $R_{6}$ in Fig. 3. Here
$R_{1}$, $R_{2}$ and $R_{3}$ are candidates for $m_{1}$. How can $m_{1}$ know
the three possible candidates among the six relay stations in the boundaries?
In new generation cellular networks, each device is equipped with GPS and
hence it knows its location. We assume that the initiator $m_{1}$ knows the
location of the receiver $m_{2}$. The assumption is very common in geographic
routing[18][31][24]. So $m_{1}$ finds out the three possible candidates in
this way: It connects $m_{1}$ and $m_{2}$. The line segment intersects the
hexagonal area at $p$. Then the relay station located at the edge where $p$ is
located and the two relay stations located at its adjacent edges are
considered as good relay candidates. In Fig. 3, $p$ is located at segment
$bc$, so $R_{2}$ located at $bc$ is considered as one of the candidates. The
two adjacent edges of $bc$ are $ab$ and $cd$, so $R_{1}$ and $R_{3}$ are also
considered as $m_{1}$’s candidates.

Figure 3: To figure out three good relay stations.

In our strategy, when a cell phone(initiator) intends to talk to a receiver in
a different cell, it calls the procedure Detection() to figure out the three
possible relay stations by its location and the destination’s location. In the
procedure, it sends out a detection message to the three relay stations. Upon
receiving the detection message, each RS calls Prediction() to calculate its
potential energy and inform the energy capacity to the initiator. The
initiator calls Selection() to find the relay station with possible highest
potential energy, hands over to it and sends its payload to the RS.

Detection()
Figure out the three candidate relay stations by its location and
destination’s location
Send detection message to the three relay stations

Figure 4: Algorithm for initiator(MU) to send detection message

Prediction()
Calculate its potential energy by formula (7)
Send the energy quantity back to the initiator

Figure 5: Algorithm of computation for potential energy

Selection()
Compare the three energy quantities from the three relay stations
Select the RS with highest energy values
Hands over to the RS
Send payload to the RS

Figure 6: Algorithm of hand over

##  IV Performance Evaluation

The performance evaluation was conducted in a simulated noiseless radio
network environment using MATLAB. We create a topology that consists of 20
hexagonal cells. A base station is located in the center of each cell and a
relay station is located at each edge. We randomly distribute a number of
mobile users to the cells. We performed a sequence of experiments in which the
number of mobile users is changed from 200 to 800 with increment of 50. For
each number of mobile users, we measure the aggregate throughput 10 times and
present the average. We also measure the average life time of the relay
stations in our experiments. The life time of a relay station is important. If
one RS runs out of energy and turns down, there will be a hole[14] and the
whole throughput will be decreased dramatically.

Our approach considers energy balance. We compare our approach(“EB by MU”)
with two other approaches. One is the scheme which considers energy balance
but the balance is determined by base station(“EB by BS”). The second one is
without considering energy balance(“No EB”).

Fig. 7 shows our evaluation of aggregate throughput. Among the three
approaches, the scheme without energy balance(“No EB”) has the lowest
throughput. This indicates that energy balance can improve aggregate
throughput because the traffic are allocated more reasonable to the relay
stations. The throughput of our method is the highest one.

Figure 7: Aggregate throughput Figure 8: Life time of relay stations

Fig. 8 shows the life time of the three schemes. The one without energy
balance(“No EB”) results in the lowest average lift time. That is because when
a relay station has busy traffic load, it runs of out energy quickly. And then
there is a void or hole. The regular solution for hole problem is to find the
landmark relay station [14]. However, the landmark has very heavy burden and
it is exhausted quickly. This makes the network much worse. Our approach has
the longest lift time because the communication between two cell phones is
always relayed by the RS with the highest potential energy. If a RS is busy
and its energy goes down significantly, the new communication will not rely on
it. Hence the balanced energy utilization makes the life time longer.

##  V Conclusion

In this paper, we present a statistical based mechanism to balance energy
utilization in cellular network. In the approach, each relay station only
needs to remain two factor to store its historical energy trend and then it is
able to predict its future energy. The message initiator selects the relay
station with potential highest energy to relay data transmission. Simulations
demonstrate that our approach results in higher throughput in the network and
longer average life time of relay stations over related approaches.

## References

  * [1]
  * [2] G. Fettweis and Z. Ernesto, _ICT energy consumption-trends and challenges_ ,Proceedings of the 11th International Symposium on Wireless Personal Multimedia Communications. Vol. 2. No. 4. 2008.
  * [3] H. Ziaul, H. Boostanimehr and V.K. Bhargava, _Green cellular networks: A survey, some research issues and challenges_ ,Communications Surveys & Tutorials, IEEE 13.4 (2011): 524-540.
  * [4] M. Gruber, et al., _EARTH: energy aware radio and network technologies_ ,ersonal, Indoor and Mobile Radio Communications, 2009 IEEE 20th International Symposium on. IEEE, 2009.
  * [5] Y. Chen, et al., _Fundamental trade-offs on green wireless networks_ ,Communications Magazine, IEEE 49.6 (2011): 30-37.
  * [6] X. wang, et al., _A survey of green mobile networks: Opportunities and challenges_ ,obile Networks and Applications 17.1 (2012): 4-20.
  * [7] M. Etoh, T. Ohya and Y. Nakayama, _Energy consumption issues on mobile network systems_ ,Applications and the Internet, 2008. SAINT 2008. International Symposium on. IEEE, 2008.
  * [8] J. Yang and Z. Fei, _Bipartite Graph Based Dynamic Spectrum Allocation for Wireless Mesh Networks_ , Distributed Computing Systems Workshops, 2008. ICDCS’08. 28th International Conference on (pp. 96-101). IEEE, Beijing, 2008.
  * [9] G. Josh, _On Relay-assisted Cellular Networks_ , thesis, Department of Electrical Engineering, Indian Institute of Technology-Mombay, June, 2010
  * [10] S. Patra, et al., _Improved genetic algorithm for channel allocation with channel borrowing in mobile computing_ ,Mobile Computing, IEEE Transactions on 5.7 (2006): 884-892.
  * [11] S. Hanly, _An algorithm for combined cell-site selection and power control to maximize cellular spread spectrum capacity_ ,Selected Areas in Communications, IEEE Journal on 13.7 (1995): 1332-1340.
  * [12] M. Soleimanipour, W. Zhuang and G.H. Freeman, _Optimal resource management in wireless multimedia wideband CDMA systems_ ,SMobile Computing, IEEE Transactions on 1.2 (2002): 143-160.
  * [13] E. Yanmaz, and O. Tonguz, _Dynamic load balancing and sharing performance of integrated wireless networks_ , Selected Areas in Communications, IEEE Journal on 22.5 (2004): 862-872.
  * [14] J. Yang, and Z. Fei, _ITGR: Intermediate Target Based Geographic Routing_ , Computer Communications and Networks (ICCCN), 2010 Proceedings of 19th International Conference on (pp. 1-6). IEEE, Zurich, Switzerland, August 2010.
  * [15] S. Alam, et al., _Energy efficient relay-assisted cellular network model using base station Switching_ , Globecom Workshops (GC Wkshps), 2012 IEEE. IEEE, 2012.
  * [16] K. XU and M. Zhou, _Energy balanced chain in IEEE 802.15. 4 low rate WPAN_ , Computing, Networking and Communications (ICNC), 2013 International Conference on. IEEE, 2013.
  * [17] K. XU and I. Howitt, _Realistic energy model based energy balanced optimization for low rate WPAN network_ , Southeastcon, 2009. SOUTHEASTCON’09. IEEE. IEEE, 2009.
  * [18] J. Yang and Z. Fei, _HDAR: Hole detection and adaptive geographic routing for ad hoc networks_ , Computer Communications and Networks (ICCCN), Proceedings of 19th International Conference. IEEE, 2010.
  * [19] G. Joshi, _On Relay-assisted Cellular Networks_ , Graduate thesis, India Institute of Technology-Bombay, 2010
  * [20] Y. Su, Y. Wang, G. Agrawal and R. Kettimuthu, _SDQuery DSI: integrating data management support with a wide area data transfer protoco_ ,Proceedings of SC13: International Conference for High Performance Computing, Networking, Storage and Analysis, 2013.
  * [21] H. Zhang, Z.Zhang and H. Dai, _On the Capacity Region of Cognitive Multiple Access over White Space Channels_ , Selected Areas in Communications, IEEE Journal on, Vol. 31, Number 11, Pages 2517–2527, (2013).
  * [22] J. Shen, P. Su and S. Cheung, _Virtual Mirror Rendering with Stationary RGB-D Cameras and Stored 3D Background_ , IEEE Transactions on Image Processing, vol. 22, issue 9, pp. 1-16.
  * [23] K. Xu, D. Tipper, P. Krishnamurthy and Y. Qian, _An Efficient Hybrid Model and Dynamic Performance Analysis for Multihop Wireless Networks_ , in Proceedings of 2013 IEEE International Conference on Computing, Networking and Communications (ICNC’13), San Diego, CA, USA, January 2013.
  * [24] J. Yang and Z. Fei, _Broadcasting with Prediction and Selective Forwarding in Vehicular Networks_ , International Journal of Distributed Sensor Networks, 2013 (2013).
  * [25] L. Li, _Evaluating Information Retrieval Systems With Multiple Non-Expert Assessors_ , (2013).
  * [26] J. Shen and S. Cheung _Layer Depth Denoising and Completion for Structured-Light RGB-D Cameras_ , IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2013), Portland, USA 2013.
  * [27] J. Shen and W. Tan, _Image-based indoor place-finder using image to plane matching_ , IEEE International Conference on Multimedia and Expo (ICME 2013), San Jose, USA 2013.
  * [28] J. Shen et al., _Automatic lip-synchronized video-self-modeling intervention for voice disorders_ , in e-Health Networking, Applications and Services (Healthcom), 2012 IEEE 14th International Conference on, 2012, pp. 244-249.
  * [29] Q. Sun, W. Yu, Q. Hao and F. Hu, _A Multi-agent Based Intelligent Sensor and Actuator Network Design for Smart House and Home Automation_ ,Journal of Sensor and Actuator Network. vol. 2, no. 3, pp. 557-588, 2013.
  * [30] J. Shen, A. Raghunathan, S. Cheung and R. Patel, _Automatic Content Generation for Video Self Modeling_ , IEEE International Conference on Multimedia and Expo (ICME 2011), Barcelona, Spain 2011.
  * [31] J. Yang and Z. Fei, _Statistical Filtering Based Broadcast Protocol for Vehicular Networks_ , 20th International Conference on Computer Communication Networks, Maui, Hawaii, USA, 2011.
  * [32] H. Zhang, Z. Zhang and H. Dai, _Gossip-based Information Spreading in Mobile Networks_ , (2012).
  * [33] H. Zhang, Z. Zhang and H. Dai, _Mobile Conductance in Sparse Networks and Mobility-Connectivity Tradeoff_ , arXiv preprint arXiv:1403.7691 (2014).
  * [34] Q. Sun, F. Hu and Q. Hao, _Mobile Target Scenario Recognition via Low-cost Pyroelectric Sensing System: Towards a Context-Enhancement Accurate Identification_ ,IEEE Transactions on Systems, Man and Cybernetics: Systems. pp. 2168-2216, 2013.
  * [35] Q. Sun, R. Ma, F. Hu and Q. Hao, _Space Encoding Based Human Activity Modeling and Situation Perception_ ,In Proceeding of 2013 IEEE International Multi-Disciplinary Conference on Cognitive Methods in Situation Awareness and Decision Support, pp. 186-189, 2013.
  * [36] Y. Wang, W. Jiang and G. Agrawal, _SciMATE: A Novel MapReduce-Like Framework for Multiple Scientific Data Formats_ ,Cluster, Cloud and Grid Computing (CCGrid), 2012 12th IEEE/ACM International Symposium on, PP 443–450, 2012.
  * [37] Y. Wang, Y. Su and G. Agrawal, _Supporting a Light-Weight Data Management Layer over HDF5_ ,Cluster, Cloud and Grid Computing (CCGrid), 2013 13th IEEE/ACM International Symposium on, PP 335-342, 2013.
  * [38] L. Li and M. Smucker, _Tolerance of Effectiveness Measures to Relevance Judging Errors_ , European conference on Information Retrieval, 2014.
  * [39] Y. Zhou, L. Li and H. Zhang, _Adaptive Learning of Region-based pLSA Model for Total Scene Annotation_ , ArXiv preprint arXiv:1311.5590 (2013).
  * [40] Y. Zhou, L. Li, T. Zhao and H. Zhang, _Region-based high-level semantics extraction with CEDD_ , Network Infrastructure and Digital Content, 2010 2nd IEEE International Conference on. IEEE, 2010.


# Binding prediction of multi-domain cellulases with a dual-CNN

Kay S. Schaller<EMAIL_ADDRESS>[  Jeppe Kari  [  Kim Borch  [  Günther H.J.
Peters  [  Peter Westh  [

###### Abstract

Cellulases hold great promise for the production of biofuels and biochemicals.
However, they are modular enzymes acting on a complex heterogeneous substrate.
Because of this complexity, the computational prediction of their catalytic
properties remains scarce, which restricts both enzyme discovery and enzyme
design. Here, we present a dual-input convolutional neural network to predict
the binding of multi-domain enzymes. This regression model outperformed
previous molecular dynamics-based methods for binding prediction for
cellulases in a fraction of the time. Also, we show that when changed to a
classification problem, the same network can be back-propagated to suggest
mutations to improve enzyme binding. A similar approach could increase our
understanding of the structure-activity relationship of enzymes, and suggest
new promising mutations for enzyme design using explainable artificial
intelligence.

###### keywords:

Cellulases, computational biology, computer modeling,protein design,protein
engineering

DTU Bioengineering] Department of Biotechnology and Biomedicine, Technical
University of Denmark [DTU Kemi] Department of Chemistry, Technical University
of Denmark DTU Bioengineering] Department of Biotechnology and Biomedicine,
Technical University of Denmark Novozymes] Novozymes A/S DTU Kemi] Department
of Chemistry, Technical University of Denmark DTU Bioengineering] Department
of Biotechnology and Biomedicine, Technical University of Denmark CAZy:
Carbohydrate-Active Enzyme, CBH: Cellobiohydrolase, CBM: Cellulose-Binding
Module, CD: Catalytic Domain, Cel: Cellulase, CNN: Convolutional Neural
Network, EG: Endoglycanase, Grad-CAM: Gradient-weighted Class Activation
Mapping GH: Glycoside Hydrolase, LIE: Linear Interaction Energy (method), MAE:
Mean-Absolute Error, MD: Molecular Dynamics, ML: Machine Learning, MSE: Mean-
Squared Error, RMSE: Root-Mean Squared Error, SGD: Stochastic Gradient
Descent, Tr: Trichoderma reesei, XAI: Explainable Artificial Intelligence,

##  1 Introduction

Cellulose is the most abundant organic compound on Earth and its utilization
promises a way towards CO2-neutral fuels, materials, and chemicals.1, 2
However, it is very recalcitrant and hard to utilize. Over millennia, Nature
developed a toolbox of enzymes to make use of this polymer. Primarily certain
fungi and bacteria are known for their biomass-degrading capabilities. For
simplicity, we will focus on the well studied and industrial very relevant
fungal cellulases.2

Figure 1: Illustration of the cellulase Cel7A from the fungus Trichoderma
reesei (TrCel7A) bound on a cellulose fibril highlighting the typical domains
of modular cellulases (left). The catalytic domain (CD) is bound to the
carbohydrate-binding module (CBM) by a flexible linker.2 Previous work showed
that the binding of the whole complex could be approximated by the binding of
the CD towards a threaded chain and binding of the CBM on the crystalline
surface (right).3, 4, 5

Cellulases are typically modular, consisting of a catalytic domain (CD) and a
carbohydrate-binding module (CBM) connected by a flexible linker (see Fig. 1).
The interactions between the enzymes and their substrate are very complex and
non-trivial to set up for simulation. This complexity hinders computational
exploration of these important enzymes compared to enzymes with a simpler mode
of interaction. One important parameter to guide enzyme engineering is their
binding affinity. However, computing of binding affinity in silico is
cumbersome and computationally expensive.

The enzyme-substrate complex of cellulases is mainly formed through
interactions of the CD towards the threaded substrate chain and the CBM on the
crystal surface. Earlier work showed that it can be approximated by the
individual domains’ binding to the substrate (see Fig. 1).3, 4, 5 This
simplifies the setup drastically, as crystal structures of the CD in complex
with threaded substrate exists and the interactions of the CBM towards the
crystal surface are well studied.2

Building up on this approximation, we recently proposed a cheap computational
method based on molecular dynamics (MD) simulation.4 While this method already
has proven useful to screen through smaller enzyme families to discover new
promising enzymes, it is still not feasible to compute properties of larger
families or explore larger sequence space.5, 6 Here, machine learning (ML)
techniques promise to capture the complexity of the underlying interactions,
providing reasonable binding estimates from either sequence or structure in an
even shorter time.7, 8, 9 Additionally, when at least partially based on real
experimental results, machine learning could not only lead to faster
predictions but also to more accurate ones than compared to the semi-empirical
MD-based 2D-LIE screening method.10, 4

Cellulases can be found in different glycoside hydrolase (GH) families and
have therefore different folds of the catalytic domain. While bacterial
cellulases can possess CBMs from different CBM families, the here investigated
fungal enzymes only have CBMs from a single family (CBM1).11, 2 These
different folds of the CDs imply that sequence and their related structural
parts are not aligned across the families. However, as they act on the same
substrate and consist of the same building blocks, their structure-function
relationship should be common. This is also indicated by the capability of the
semi-empirical 2D-LIE to predict across different families by utilizing common
scaling parameters.4 Because of this, we focus on structure-based machine
learning, which hopefully allows the model to easier transfer knowledge across
families. Additionally, there are plenty of crystal structures available for
the well-studied cellulases, which makes modeling of their structure very
straightforward, even before the advent of Alphafold (2.0).4, 12

Within machine learning, there is a wide array of methods, but in this work,
because we want to start from structure, methods originating from image
recognition are well suited for the task. Especially, convolutional neural
networks (CNN) have been shown to be capable of capturing underlying patterns
in structural data, such as medical imaging data.13 Earlier work also showed
that it is possible to generalize across different protein structures and
predict protein-ligand binding via CNN.14, 15, 16 However, those efforts
focused on globular enzymes binding to small molecules. As cellulases are
multi-domain enzymes binding to an interfacial substrate, those approaches are
not readily transferable.

##  2 Experimental procedures

###  2.1 Data set & data curation

In two recent studies, we performed virtual screenings with MD-based binding
estimates for two cellulase families, GH6 and GH7.5, 6 Family GH7 consisted of
approx. 400 representative enzymes, family GH6 of about 1’200. Additionally,
we earlier published a (relatively) large dataset of ca. 100 cellulases with
experimental binding values.3 For a more efficient training, all binding
energies were linearly scaled with a minimum value of $-13$ kJ/mol and a
maximum value of $17$ kJ/mol, resulting in values well between $0$-$1$ for all
datasets. The structures were aligned per domain type and split into protein
and organic molecules. Moleculekit was used to voxelize the structures
according to earlier work by Jimenez and co-workers with a resolution of
$\mathrm{1}$ Å.14 The original publication used 8 different channels. For the
protein parts (CD/CBM), we dropped the ”metal ion”-encoding channels as all of
our homology-modeled structures are free of ions, thus resulting in 7 channels
for the protein part of each domain (see Fig. 2). For the organic molecule
part (polymeric ligand/crystal surface), we additionally dropped the channels
encoding aromatic, hydrophobic, positive ionizable, and negative ionizable
atoms, because cellulose does not have atoms with those attributes. This led
to 3 channels for the cellulose part of each domain. In total, this gave us 10
channels per domain, 7 from the protein part (CD/CBM) and 3 from the cellulase
part (bound ligand/crystal surface, see Fig. 2).

Figure 2: Input example for the enzyme TrCel7A. 10 channels were employed for
both domains (CD with polymeric ligand and CBM on surface). 7 channels encode
the protein (CD/CBM) and 3 encoding the organic molecules (polymeric
ligand/surface). The CBM of TrCel7A does not have negative ionizable atoms and
the 6th input channel of that domain is therefore empty.
Figure 3: Architecture of the dual-input CNN for binding prediction of multi-
domain cellulases. For regression a single output node is used, while for
classification 2 nodes were used to distinguish between ”good” and ”weak”
binders.

As inputs, a $\mathrm{50}$ Å cube from the voxelized CD data and a
$\mathrm{30}$ Å from the CBM data was taken. In all cases, an $80\%/20\%$
train-test split was used. Enzyme families were proportionally stratified
across the test and training set. To increase generalization and allow
positional invariance, the training data was augmented. Random translations of
the window taken as input cube in all three dimensions (up to 12 Å in all
directions) were done. Additionally, random rotations for all possible 24 cube
orientations, and random flips in all three dimensions were performed for the
CD and CBM input. Accounting for permutations between the CD and CBM input,
this results in millions of possible combinations. On-the-fly augmentation was
implemented and tested but led to poor computational performance. Therefore,
all datasets were split beforehand, and a fixed number of augmentation were
done on the training data. For training on the larger dataset based on GH6 and
GH7, we generated $5$ augmented inputs per data point; for the smaller dataset
based on the experimental dataset, we generated $20$ augmented inputs per data
point. Ideally speaking, even the larger simulated set is still regarded as
small for deep learning applications, but biological data is costly to obtain.

###  2.2 Network

The applied CNN architecture was inspired by literature.14, 15 Those
literature models, however, only treat globular (single domain) enzymes. To be
able to predict binding for our two-domain enzymes such as cellulases, we
employ two input and two convolutional parts, one for each domain (see Fig.
3). The CD input, starting from a $\mathrm{50}$ Å, was fed into three
3D-convolutional layers (64, 128, and 256 filters, respectively). The smaller
CBM input was convoluted with two layers (64 and 128 filters, respectively).
For all convolutional layers, a kernel size of 3 was used. All convolutional
layers were followed by maximum pooling layers. The convolutional parts were
followed by a global maximum 3D-pooling layer, and a concatenation between CD
and CBM feeds. This was followed by dense layers with 1000, 496, and 200
neurons (multiples of 8 for more efficient GPU training) with $\mathrm{ReLu}$
activation. For regression training, this was fed into a single output neuron
with a $\mathrm{sigmoid}$ activation function. For classification, this was
swapped with an output layer of $2$ neurons with $\mathrm{SoftMax}$
activation.

###  2.3 Training

Because we have datasets with different accuracy (simulated and experimental
target data), we employed fine-tuning. Initially, we learned on the large
dataset based on simulated data, and after that refined the model while
learning on the more accurate but smaller set with experimental target values.
Later, we used this pre-trained regression model and transferred it to a
classification model for class-specific guided back-propagation.

#### Initial Learning on Simulated Values

Different hyperparameters were investigated. The final learning used
stochastic gradient descent (SGD) with no momentum and a fixed learning rate
of $10^{-4}$ for optimization.17 Mean-squared error (MSE) was used for loss,
and mean-absolute error (MAE) was tracked for performance checks. A batch size
of $256$ with pre-fetching and interleaved file reading was used for all
trainings. Batch normalization was applied to the convolutional layers, and a
dropout of $0.6$ was used in the dense layers.

#### Fine-Tuning on Experimental Values

Only a few things were changed compared to the initial learning: The learning
rate was decreased to $10^{-5}$, and the dense layers were kept frozen during
training.

#### Transfer-Learning for Classification Model

The output layer was changed from a single regression neuron to $2$ neurons
for classification. An arbitrary energy limit was introduced at $0$ kJ/mol to
split the dataset into a ”good binder” and a ”weak binder” class. Continuous
target values were on-the-fly converted to sparse encoded vectors. Categorical
cross-entropy loss and accuracy were used to guide the learning. The loss was
class-weighted to account for data imbalance in the different classes. A
slower learning rate of $10^{-6}$ was used while all dense layers were kept
frozen.

#### Engineering Guided by Grad-CAM

The trained classification model was taken, and gradient-weighted class
activation mapping (Grad-CAM) was performed to investigate the activation of
the input layers per class.18 TrCel7A-$\mathrm{W38A}$ and
TrCel7A-$\mathrm{W40A}$, two known mutants of TrCel7A were taken as a show-
case. It is known that the residues $\mathrm{W38}$ and $\mathrm{W40}$ are
critical for binding and that their mutation to alanine leads to a drastic
decrease in binding energy.3, 19

#### Software

The deep learning model was built with Tensorflow (2.4) using python (3.8.11),
CUDA (11.3), scikit-learn (0.24.2).20, 21 The module 3D-GuidedGradCAM was used
to perform guided back-propagation.18

##  3 Results

Figure 4: Overview over different learning stages with the used datasets,
exemplary learning curves, and prediction plots for the validation set. In the
prediction plots (right), the performance of the previous 2D-LIE method
($\pm$RMSE) is shown as dashed lines.4 Top: Initial learning of the regression
model on the larger dataset with simulated values. Center: Fine-tuning of
regression model on the smaller dataset with experimental values. Bottom:
Transfer learning of the classification model based on the pre-trained
regression model. Guided Grad-CAM can be performed, and subsequent class
discriminant analysis can elucidate mutations for improved binding.

#### Initial Learning on Simulated Values

Learning on the larger dataset based on simulated target values was performed
over 1000 epochs (see top of Fig. 4). No increase in validation performance
was observed at the end, suggesting that no overfitting has yet been reached.
The final model resulted in a root-mean-squared error (RMSE) of $2.09$ kJ/mol,
only slightly higher than the estimated $1.91$ kJ/mol for the initial
simulation method at longer simulation times ($10$ ns instead of the $1$ ns,
that was used for the generation of the dataset), which could be close to the
accuracy limit of the dataset.

#### Fine Tuning on Experimental Values

Learning was performed over 500 epochs (see center of Fig. 4). Again, no
increase in validation performance was observed at the end, suggesting that
overfitting has not been reached yet. The final model resulted in an RMSE of
$1.53$ kJ/mol, outperforming the previous method 2D-LIE (RMSE of $1.98$
kJ/mol). This performance is comparable to related methods for globular
enzymes from literature, despite predicting binding of more complex modular
enzymes. 15, 14, 16, 9 The regression model seems to work well for all
included families (GH5, GH6, GH7, and GH45) and irrespective of modularity
($\pm$CBM).

#### Transfer-Learning for Classification Model

Learning was performed over 100 epochs leading to a perfect categorical
accuracy on the validation set (see bottom of Fig. 4). The jumps in the
accuracy suggest that using a different, more stable optimizer for the
training could be advantageous.

Figure 5: Class discriminant analysis of two TrCel7A variants (W38A and W40A)
with Grad-CAM. Position 38 is highlighted in orange, position 40 in gold. Both
residues are found in the tunnel entrance (right) and are known to stack with
the threaded cellulose chain. Inlets (left) show ”good binder” class
activation (top 25%) as green isomeshes for the wildtype TrCel7A and the two
variants. The presence of the tryptophan residue leads to class-specific
activation, as seen by the appearance of the isomeshes.

#### Engineering Guided by Grad-CAM

The wildtype (TrCel7A) is found in the ”good binder” class ($<0$ kJ/mol), both
according to the ground truth and predicted by the classification model. The
point mutations are not impactful enough to change the class, however, the
change in activation can still be interpreted. TrCel7A,
TrCel7A-$\mathrm{W38A}$, and TrCel7A-$\mathrm{W40A}$ were inferred onto the
model, and Grad-CAM was used to perform class discriminant analysis of the
activation (see Fig. 5). The presence of the tryptophans leads to an increased
activation for the ”good binder” class. In other words, the model becomes more
sure that the enzyme binds well.

##  4 Discussion

The structure-activity relationship of proteins and their substrate is non-
linear and very complex in general. Thus, they pose an ideal target for
machine learning if a feasible encoding for the problem can be found. One
possible approach is the voxelization of the structure. Ideally, this readily
allows to capture fold-independent properties as it is solely structure and
not sequence-based.14, 7 In this work, we presented machine-learning models
for multi-domain enzymes in general and cellulases in particular based on
voxelized structures as inputs. Cellulases are an important class of enzymes
relevant for the conversion of biomass to fuels and chemicals. Their redesign
is relevant for improved activity under industrial conditions, and one
important descriptor for their engineering is their binding affinity. However,
the computational assessment of their binding is non-trivial as their enzyme-
substrate interaction is very complex. Therefore, we developed a method to
predict their binding affinity to guide future enzyme engineering.

We used simulated binding data from two previous studies of family GH6 and GH7
to pre-train a regression model.6, 5 Subsequently, we used a smaller set of
experimental values4 to fine-tune the regression model, resulting in a model
with an estimated RMSE of $1.53$ kJ/mol, outperforming previous computational
screening methods for cellulases. Once trained, this regression model is also
orders of magnitude faster than the previous MD-based method. The trained
regression models were transferred to a classification model to allow class-
specific guided back-propagation. The classification model itself is not that
practical, but it allows to use of explainable artificial intelligence (XAI)
methods developed for classification problems.22, 23, 24 Two mutant enzymes
with poor binding were class-specific back-propagated, and the found densities
suggested reintroducing the mutationally removed aromatic residues. While this
reversal to increase binding once again is trivial, the same back-propagation
approach could be used to find other unknown mutational suggestions to improve
binding. Future work could directly use networks with simultaneous regression
and classification output nodes to avoid retraining and additional hyper-
parameter search. More general, regression and classification models like
these could be investigated with additional XAI methods to learn more about
the structure-activity relationship of enzymes. Approaches like these could
make it feasible to search vast sequence space for interfacial enzymes such as
cellulases and therefore push enzyme discovery and design to new limits. Using
XAI methods could allow us to understand more of the hidden structure-activity
relationship within the ML models and lead to a new understanding of enzyme-
substrate interaction.

###  4.1 Data availability

All input data is based on public datasets as described in the referenced
papers. All employed software is open-source and all necessary information is
listed.

###  4.2 Conflict of interest

K.B. works for Novozymes A/S, a major enzyme-producing company.

###  4.3 Acknowledgments

The simulations were carried out at the high-performance cluster at the
Technical University of Denmark. This work was supported by the Independent
Research Fund Denmark (Grant No. 8022-00165B) and the Novo Nordisk Foundation
(Grant Nos. NNF15OC0016606 and NNF17SA0028392).

## References

  * Ragauskas et al. 2006 Ragauskas, A. J.; Williams, C. K.; Davison, B. H.; Britovsek, G.; Cairney, J.; Eckert, C. A.; Frederick, W. J. J.; Hallett, J. P.; Leak, D. J.; Liotta, C. L.; Mielenz, J. R.; Murphy, R.; Templer, R.; Tschaplinski, T. The Path Forward for Biofuels and Biomaterials. _Science_ 2006, _311_ , 484–489
  * Payne et al. 2015 Payne, C. M.; Knott, B. C.; Mayes, H. B.; Hansson, H.; Himmel, M. E.; Sandgren, M.; Ståhlberg, J.; Beckham, G. T. Fungal Cellulases. _Chemical Reviews_ 2015, _115_ , 1308–1448
  * Kari et al. 2021 Kari, J. et al. Physical Constraints and Functional Plasticity of Cellulases. _Nature Communications_ 2021, _12_ , 3847
  * Schaller et al. 2021 Schaller, K. S.; Kari, J.; Molina, G. A.; Tidemand, K. D.; Borch, K.; Peters, G. H. J.; Westh, P. Computing Cellulase Kinetics with a Two-Domain Linear Interaction Energy Approach. _ACS Omega_ 2021, _6_ , 1547–1555
  * Schaller et al. 2022 Schaller, K. S.; Molina, G. A.; Kari, J.; Schiano-di Cola, C.; Sørensen, T. H.; Borch, K.; Peters, G. H.; Westh, P. Virtual Bioprospecting of Interfacial Enzymes: Relating Sequence and Kinetics. _ACS Catal._ 2022, 7427–7435
  * Kay Sven Schaller et al. 2022 Kay Sven Schaller,; Jeppe Kari,; Gustavo Avelar Molina,; Kim Borch,; Günther H.J. Peter,; Peter Westh, Virtual Screening of GH Family 6\. _Manuscript in preparation_ 2022, _0_ , 1–2
  * Greener et al. 2021 Greener, J. G.; Kandathil, S. M.; Moffat, L.; Jones, D. T. A Guide to Machine Learning for Biologists. _Nature Reviews Molecular Cell Biology 2021_ 2021, 1–16
  * Mazurenko et al. 2019 Mazurenko, S.; Prokop, Z.; Damborsky, J. Machine Learning in Enzyme Engineering. _ACS Catalysis_ 2019, _10_ , 1210–1223
  * Ellingson et al. 2020 Ellingson, S. R.; Davis, B.; Allen, J. Machine Learning and Ligand Binding predictions: A Review of Data, Methods, and Obstacles. _Biochimica et Biophysica Acta - General Subjects_ 2020, _1864_ , 129545
  * Noé et al. 2020 Noé, F.; Tkatchenko, A.; Müller, K.-R.; Clementi, C. Machine Learning for Molecular Simulation. _Annual Review of Physical Chemistry_ 2020, _71_ , 361–390
  * Henrissat 1991 Henrissat, B. A Classification of Glycosyl Hydrolases Based on Amino Acid Sequence Similarities. _Biochemical Journal_ 1991, _280_ , 309–316
  * Jumper et al. 2021 Jumper, J. et al. Highly accurate protein structure prediction with AlphaFold. _Nature 2021_ 2021, 1–7
  * Yamashita et al. 2018 Yamashita, R.; Nishio, M.; Do, R. K. G.; Togashi, K. Convolutional neural networks: an overview and application in radiology. _Insights into Imaging 2018 9:4_ 2018, _9_ , 611–629
  * Jiménez et al. 2018 Jiménez, J.; Škalič, M.; Martínez-Rosell, G.; De Fabritiis, G. KDEEP: Protein-Ligand Absolute Binding Affinity Prediction via 3D-Convolutional Neural Networks. _Journal of Chemical Information and Modeling_ 2018, _58_ , 287–296
  * Stepniewska-Dziubinska et al. 2018 Stepniewska-Dziubinska, M. M.; Zielenkiewicz, P.; Siedlecki, P. Development and Evaluation of a Deep Learning Model for Protein–Ligand Binding Affinity Prediction. _Bioinformatics_ 2018, _34_ , 3666–3674
  * Jones et al. 2021 Jones, D.; Kim, H.; Zhang, X.; Zemla, A.; Stevenson, G.; Bennett, W. F. D.; Kirshner, D.; Wong, S. E.; Lightstone, F. C.; Allen, J. E. Improved Protein–Ligand Binding Affinity Prediction with Structure-Based Deep Fusion Inference. _Journal of Chemical Information and Modeling_ 2021, _61_ , 1583–1592
  * Bottou and Bousquet 2012 Bottou, L.; Bousquet, O. _Optimization for Machine Learning_ ; MIT Press, 2012; Vol. 2; p 152
  * Selvaraju et al. 2017 Selvaraju, R. R.; Cogswell, M.; Das, A.; Vedantam, R.; Parikh, D.; Batra, D. Grad-CAM: Visual Explanations from Deep Networks via Gradient-Based Localization. _Proceedings of the IEEE International Conference on Computer Vision_ 2017, _2017_ , 618–626
  * Røjel et al. 2020 Røjel, N.; Kari, J.; Sørensen, T. H.; Badino, S. F.; Morth, J. P.; Schaller, K.; Cavaleiro, A. M.; Borch, K.; Westh, P. Substrate Binding in the Processive Cellulase Cel7A: Transition State of Complexation and Roles of Conserved Tryptophan Residues. _Journal of Biological Chemistry_ 2020, _295_ , 1454–1463
  * Martín Abadi et al. 2015 Martín Abadi, et al. TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems. 2015
  * Pedregosa et al. 2011 Pedregosa, F. et al. Scikit-learn: Machine Learning in Python. 2011; http://scikit-learn.sourceforge.net.
  * Emmert-Streib et al. 2020 Emmert-Streib, F.; Yli-Harja, O.; Dehmer, M. Explainable Artificial Intelligence and Machine Learning: A Reality Rooted Perspective. _WIREs Data Mining and Knowledge Discovery_ 2020, _10_ , e1368
  * Angelov et al. 2021 Angelov, P. P.; Soares, E. A.; Jiang, R.; Arnold, N. I.; Atkinson, P. M. Explainable Artificial Intelligence: An Analytical Review. _WIREs Data Mining and Knowledge Discovery_ 2021, _11_ , e1424
  * Confalonieri et al. 2021 Confalonieri, R.; Coba, L.; Wagner, B.; Besold, T. R. A Historical Perspective of Explainable Artificial Intelligence. _WIREs Data Mining and Knowledge Discovery_ 2021, _11_ , e1391